Strong increase in thawing of subsea permafrost in the 22nd century caused by anthropogenic climate change

Most earth system models (ESMs) neglect climate feedbacks arising from carbon release from thawing permafrost, especially from thawing of subsea permafrost (SSPF). To assess the fate of SSPF in the next 1000 years, we implemented SSPF into JSBACH, the land component of the Max Planck Institute Earth System Model (MPI-ESM). This is the first implementation of SSPF processes in an ESM component. We investigate three extended scenarios from the 6th phase of the Coupled Model Intercomparison Project (CMIP6). In the 21st century only small differences are found among the scenarios, but in the upper-end emission scenario SSP5–8.5 (shared socio-economic pathway), especially in the 22nd century, SSPF ice melting is more than 15 times faster than in the pre-industrial period. In this scenario about 35 % of total SSPF volume and 34 % of SSPF area are lost by the year 3000 due to climatic changes. In the more moderate scenarios, the melting rate maximally exceeds that of pre-industrial times by a factor of 4, and the climate change induced SSPF loss (volume and area) by the year 3000 does not exceed 14 %. Our results suggest that the rate of melting of SSPF ice is related to the length of the local open-water season and thus that the easily observable sea ice concentration may be used as a proxy for the change in SSPF.

1 Introduction

More than 1300 Pg of carbon are stored in the permafrost soils (perennially frozen soils) of the Arctic (Hugelius et al., 2014). During the present interglacial period, the Holocene, microbiological activity in partially thawed soils degraded a fraction of the stored organic carbon and released it to the atmosphere. Enhanced warming during the Anthropocene has in recent years accelerated this carbon release (Schuur et al., 2015). Most studies have been focused on the thawing of the terrestrial permafrost (Koven et al., 2015; Kleinen and Brovkin, 2018; Turetsky et al., 2020). Since the Last Glacial Maximum (LGM) about 3.5×10⁶ km² (Sayedi et al., 2020) of permafrost soils were submerged by the sea level rising by about 120 m. These submerged permafrost sediments (subsea permafrost, SSPF) now form the major part of the Arctic Shelf. Sayedi et al. (2020) suggest that about 500 Pg of carbon in the form of organic carbon and methane gas may be trapped in the SSPF and become available for microbial decomposition as SSPF thaws. This estimate is however highly uncertain due to the sparsity of measurements (Schuur et al., 2015). Since the benthic temperatures are above the freezing point and their inter-annual variability is small, the submergence below seawater causes a slow but continuous thawing of the submerged permafrost sediments from the top. Unless the Earth enters a new glacial state, the SSPF will thus ultimately be thawed away at some time. Anthropogenic climate change may accelerate the thawing and thus the release of carbon.

The quantitative projection of the future climate is often done using earth system models (ESMs) which include the carbon cycle. However, among the ESMs participating in the 6th phase of the Coupled Model Inter-comparison Project (CMIP6; Eyring et al., 2016), only a few include the dynamics of the terrestrial permafrost (Lawrence et al., 2019), and none of them included subsea permafrost. These models thus lack potentially important feedbacks from the carbon cycle on the global climate. Several modelling studies have aimed at projecting the fate of terrestrial permafrost (Kleinen and Brovkin, 2018; McGuire et al., 2018). A few studies hindcasted subsea permafrost regionally (Shakhova et al., 2009; Nicolsky and Shakhova, 2010; Nicolsky et al., 2012; Angelopoulos et al., 2019) or at the pan-Arctic scale (Overduin et al., 2019). Extrapolating observed temperature trends, Dmitrenko et al. (2011) constructed a future scenario for the SSPF in the Laptev Sea and on the East Siberian Shelf.

In this study we present the first steps to include subsea permafrost in the Max Planck Institute Earth System Model (MPI-ESM) to investigate the magnitude of SSPF–climate feedbacks in the future. We explore the possible range of pan-Arctic SSPF thawing by applying climate forcings from several CMIP6 scenarios to the Arctic SSPF areas.

2 Methods

The land model JSBACH (Reick et al., 2013; Brovkin et al., 2013) is the land component of the MPI-ESM (Giorgetta et al., 2013). The JSBACH version in the MPI-ESM 1.2 (Mauritsen et al., 2019) includes a multi-layer soil hydrology (Hagemann and Stacke, 2015) and the multi-compartment soil carbon model YASSO (Tuomi et al., 2009; Goll et al., 2015). The version used in this study was extended with a model for freezing soil water (Ekici et al., 2014), the dynamical inundation model TOPMODEL (Beven and Kirkby, 1979; de Vrese et al., 2021), and a methane release model (Riley et al., 2011) as described in Kleinen et al. (2020). Furthermore, de Vrese et al. (2021) added vertically layered soil carbon to study the organic soil development in permafrost regions. Still modifications were needed to be able to simulate SSPF.

2.1 Submerging

To implement SSPF processes in an ESM land model, it is necessary to consider both land and subsea points which experience very different boundary conditions and have different active processes.

For the subsea points MPI-OM, the ocean component of MPI-ESM (see, for example, Jungclaus et al., 2006, Jungclaus et al., 2013, and Mauritsen et al., 2019, for details), benthic temperatures were used as an upper boundary condition instead of atmospheric temperatures as are used for terrestrial points. For subsea points, the radiative and hydrological forcings were switched off. The sediment pore space is assumed to be constantly filled with a mixture of water and ice, and there is no advection of water.

Plant growth and phenology were disabled for subsea points, leaving these points with a thermally active and carbon-decomposing sediment model.

2.2 Geothermal heat flux

In most experiments JSBACH was forced by a geothermal heat flux from beneath. The geothermal heat flux was based on the geologically based data set presented in Davies (2013). This data set was processed for the SuPerMAP model (Overduin et al., 2019) and converted from a flux (units W m⁻²) to a geothermal temperature gradient (units K m⁻¹), assuming the average heat conductivity of the sediment/water/ice mixture of the lowest layers (i.e. low porosity) to be the same as that of the JSBACH bedrock: 2 W K⁻¹ m⁻¹. The resulting pan-Arctic geothermal temperature gradient (Fig. S1) ranges from 19 to 115 mK m⁻¹ (average 35 mK m⁻¹). A consequence of this implementation is that the temperature of the lowest model layer, layer n, cannot change arbitrarily. For each model point $T_n=T_{n-\mathrm{1}}+\frac{d_{n-\mathrm{1}}+d_n}{\mathrm{2}}$GG, where T_i and d_i are the temperature and thickness of layer i, respectively, n is the number of vertical layers, and GG is the geothermal temperature gradient. This implementation was for a single point version of JSBACH presented in Puglini et al. (2020).

2.3 Salinity

SuPerMAP incorporates salt into the newly deposited sediments during periods when the considered model points are inundated into the ocean, i.e. during the interglacials (Overduin et al., 2019). Since the model points have individual inundation histories, the sediments thus have horizontally and vertically variable salinities. The salt distribution from SuPerMAP was preserved in JSBACH in the form of a local freezing point depression: T_melt = −0.054 K kg g⁻¹ S, where S is the salinity (in units of g kg⁻¹). Diffusion of salt as modelled by, for example, Angelopoulos et al. (2019) and Malakhova and Eliseev (2020) is a very slow process (Harrison and Osterkamp, 1978, report a diffusion coefficient for salinity about 4 orders of magnitude lower than for temperature). Due to the comparatively short timescale covered in the present study, diffusion of salt was not implemented.

Table 1 Performed experiments. The forcing for the three SSPs is identical for the period 1850–2009. ^*: forcing repeated cyclicly. “Geoth”: geothermal heat flux included. “Freezing layers”: freezing of porewater allowed for these layers. “Sedi”: sediment heat capacity and conductivity from “St” standard JSBACH or “Mod” modified as described in the discussion.

Despite the slowness of the salinity diffusion, the upper shallow sediment layers may be intruded on short timescales. This causes the model artefact that the upper-layer porewater may seasonally freeze – especially during the winter season – if it is submerged under (saline) liquid sub-zero-temperature seawater. In reality this porewater would, through diffusive exchange with the overlying seawater, be sufficiently saline to prevent freezing. To prevent this buffering of energy in a non-physical refreezing of the pore water in the upper layers, the active freezing of porewater was disabled for subsea points. Additional experiments were performed to test the consequences of this hypothesis.

Implementing the freezing point depression due to salt allows two different definitions of SSPF – according to temperature or according to ice content. Here, we refer to SSPF only when the sediments are frozen; i.e. it is not sufficient to have sub-zero temperatures.

2.4 Model set-up and experiments

JSBACH was run on the T63 horizontal grid with a resolution of $\approx \mathrm{1.9}\times \mathrm{1.9}$^∘ or $\approx \mathrm{50}\times \mathrm{200}$ km at the Arctic latitudes around 75^∘ N where most of the SSPF is found. In this set-up, 407 points are located on the modern Arctic shelf at a distance of ≤20 km from a point in the SuPerMAP model (Overduin et al., 2019). In the SuPerMAP set-up, points with water with depths <150 m were selected since these potentially hold SSPF. In this JSBACH set-up a subsea point on average represents an area of about 11 400 km². Vertically, 22 layers were used to represent the upper 1000 m of sediment with layer thicknesses between 6.5 cm and 300 m. Generally the layer thickness increases exponentially with increasing depth, but additional layers were added to the depth range 0.5–3 m. This vertical set-up is a downward extension of the 18 soil layer model used by de Vrese et al. (2021).

Sediment porosity and water depth were adopted from the SuPerMAP points used. The porosity decreases exponentially with increasing depth from ≈0.4 in the upper layer to ≈0.1 in the lowest and is the same for all geographical points.

The thermal sediment properties were assumed to be those of the standard JSBACH “bedrock” for which the effective diffusion is calculated as a mixture of rock and water weighted with the porosity (Ekici et al., 2014; de Vrese et al., 2021).

Experiments were conducted following each of the three CMIP6 SSP (shared socio-economic pathway) scenarios, SSP1-2.6, SSP2-4.5, and SSP5-8.5, to evaluate the fate of SSPF ice in scenarios with low, intermediate, and upper-end future global warming, respectively.

Since SSPF is never in equilibrium with the boundary conditions, two additional control experiments were conducted, referring to a pre-industrial and a present-day climate. These two experiments are used to assess SSPF thawing in a hypothetical world without anthropogenic climate change and were cyclically forced with 24 years of forcing data from 1850 to 1873 (pre-industrial) and 1986 to 2009 (present day). Further experiments were performed to assess the consequences of some of the model assumptions.

Most of the performed experiments (Table 1) cover at least the period 1850 to 3000. The exceptions are the present-day control simulation which branches off from SSP1-2.6 in 2010 and thus covers only the years 2010 to 3000 and some of the sensitivity runs. The three main experiments, pmt_ssp126, pmt_ssp245, and pmt_ssp585, are identical for the period 1850 to 2009.

2.5 Initial conditions

The model was initialized with sediment ice content, temperature, and salinity from SuPerMAP (Overduin et al., 2019). For each point the data from the SuPerMAP point nearest to the centre of the respective JSBACH cell was used.

Figure 1 (a) Length of initial (1850) SSPF ice column. Reddish colour marks points entirely without ice in the sediments. (b) Vertical distribution of initial ice concentration in sediment pore space (SSPF ice saturation). Panel (b) shows the spatial average over the area containing SSPF ice in 1850.

We took advantage of the elaborated spin-up protocol applied to the SuPerMAP model (Overduin et al., 2019) which uses 450 kyr average temperatures from the intermediate complexity climate model CLIMBER-2 and geothermal heat flux to calculate steady-state profiles of temperature and SSPF ice for each point. This is longer than the ≈100 kyr estimated by Malakhova and Eliseev (2017) to be necessary for a proper equilibration of the deep sediments. This steady-state solution was used as initialization for 50 kyr long runs using idealized benthic temperatures (ranging from −1.8 to 0 ^∘C, dependent on water depth) and the point individual transient submergence history. To minimize artificial initial effects arising due to the transition from SuPerMAP (with its idealized upper boundary condition) to JSBACH (using the MPI-ESM data as forcing) the last 23.5 kyr (since the LGM) of the SuPerMAP run were rerun using a linear interpolation between the idealized benthic temperatures and the temporal average of the MPI-ESM CMIP6 pre-industrial benthic temperatures, defined as the period 1850–1873. For the last 1000 years, SuPerMAP was forced solely with the pre-industrial temperatures. Even though this approach largely eliminates transient effects, additionally ≈500 km³ SSPF ice melts until the year 2000 compared to what could be expected from the control experiment. More than 80 % of the additional melting takes place until the year 1900, which is also about the time when the melting in the climate-driven scenarios start exceeding that of the control simulation. Since this effect is small and our focus is on the future scenarios, we regard these minor artefacts as acceptable and without influence on our main results.

The SuPerMAP output at 2 m vertical resolution was assumed to be valid at the layer midpoint and interpolated linearly to the JSBACH layers. For JSBACH layers above 1 m, the values for the upper SuPerMAP layer were assumed representative. Since the seasonal temperature cycle penetrates further down (Fig. S13), the results are not sensitive to this initial assumption.

The ice saturations given by SuPerMAP were converted to “ice column lengths” for each layer i according to ice_i = $\frac{{\mathrm{sat}}_i}{{\mathrm{por}}_i{d}_i}$, where ice_i is the ice column length, sat_i is the ice saturation, por_i is the porosity, and d_i is the layer thickness. The initial ice column length accumulated over all layers (Fig. 1a) ranges from 0 to over 180 m. The initial area of the grid cells where SSPF ice is found at any depth in the sediments (initial SSPF area) is $\approx \mathrm{2.89}\times {\mathrm{10}}^{\mathrm{6}}$ km² (242 of the 407 points), where $\approx \mathrm{0.77}\times {\mathrm{10}}^{\mathrm{6}}$ km² is less than 10 m thick. The total initial ice volume is $\approx \mathrm{153}\times {\mathrm{10}}^{\mathrm{3}}$ km³. Vertically, the ice saturation (Fig. 1b) is highest at ≈6 m depth with a secondary peak at ≈75 m. The SSPF ice saturation is closely anti-correlated with the embedded salt, consistent with the freezing point depression and the inundation history of the last two glacial cycles.

2.6 Boundary conditions

The subsea points of JSBACH were forced with monthly mean benthic temperatures as upper boundary conditions. These were adopted from the MPI-ESM runs presented in detail in Kleinen et al. (2021) extending the CMIP6 scenarios beyond the year 3000. Until the year 2500 these runs were forced by the CO₂ concentration scenarios presented in Meinshausen et al. (2020). After the year 2500 MPI-ESM was forced by CO₂ concentrations from the intermediate complexity climate model CLIMBER-2 (Brovkin et al., 2012; Kleinen et al., 2016; Petoukhov et al., 2000). The temperature of the lowest of the 40 fixed-depth MPI-OM ocean model layers above the ocean floor was assumed to be the benthic temperature. This was interpolated from the GR30 grid (resolution varying within the SSPF area from ≈40 km near Greenland to ≈200 km on the East Siberian Arctic Shelf) onto the T63 grid using CDO remapcon (Schulzweida, 2019). Despite the rather coarse oceanic resolution in the Russian Arctic, MPI-OM reproduces the trend in the coastal Laptev Sea benthic temperature reported by Dmitrenko et al. (2011) remarkably well (2.1 K reported vs. 2.15 K in MPI-OM during the period 1985 to 2009, not shown). Our data however suggest that this strong trend is confined to the Laptev Sea area and temporally limited to a period not much longer than the reporting period of Dmitrenko et al. (2011). The forcing of the different scenarios diverge from 2010 onward (not from 2016 as described in the CMIP6 protocol, Eyring et al., 2016) for technical reasons (see Kleinen et al., 2021).

Figure 2 Benthic temperature. Lines are yearly averages over the area containing SSPF ice in 1850. Shaded areas show decadal means of the extreme months for each year and are combined RGB-wise.

Three scenarios were used: SSP1-2.6, SSP2-4.5, and SSP5-8.5. In all scenarios the benthic temperature rises from around the year 2000 and stabilizes at a new essentially constant level (Fig. 2). In SSP1-2.6 the rise only lasts a few decades and is in total <1 K; in SSP2-4.5 the temperature in total rises ≈2 K, mainly before the year 2150. For these scenarios, the largest temperature increases are found in the Barents Sea and Baffin Bay (Fig. S5), closely relating to the atmospheric temperature rise (Fig. S4), in which not much SSPF ice is found (Fig. 1a). The SSP5-8.5 benthic temperature in total rises almost 10 K with the strongest trend in the 22nd century (Fig. 2). For the individual model points the increase is between 7 and 16 K and is highest in Baffin Bay and lowest in the Laptev Sea (Fig. S5, right). In all scenarios the benthic temperatures increase less than average for points which contain SSPF, and the highest temperature rise is found in coastal waters.

Figure 3 Volume of pan-Arctic SSPF ice for the scenario and control experiments. Grey shading marks the historical period.

3 Results

The amount of SSPF ice decreases smoothly over the entire experiment period (Fig. 3) in all scenario and control experiments, however at very different rates. In the experiment forced with a pre-industrial forcing (pmt_pre) ice melts at an essentially constant rate of ≈7.5 km³ a⁻¹ over the entire period. Virtually the same behaviour at rates of 9.1 and 11.7 km³ a⁻¹ is seen in the present-day (pmt_curr) and SSP1-2.6 (pmt_ssp126) experiments, respectively. The experiments forced with the two SSPs with stronger warming, pmt_ssp245 and pmt_ssp585, result in an accelerating ice melting from the present day to 50–100 years in the future and thereafter in slowly decreasing melting rates. The average melting rate for the pmt_ssp245 is 17.7 km³ a⁻¹ and thus not too much higher than for the stabilization scenarios, pmt_curr and pmt_ssp126. However, the pmt_ssp585 has an average melting rate of 54.6 km³ a⁻¹, and thus on average more than 3 times as much ice is lost, even compared to the intermediate forcing scenario pmt_ssp245. The difference between the scenarios is minor in the 21st century, but the melting increases dramatically in the 22nd century in the pmt_ssp245 simulation and much more so in the pmt_ssp585 simulation. The ratio between the SSPF ice meltings in a scenario run and in the control run (pmt_pre) emphasizes the influence of climate change on the SSPF. During the 22nd century this ratio averages to 15.8 for pmt_ssp585 (Fig. 4), in which on average 116 km³ a⁻¹ of ice is lost. pmt_ssp245 melts 3–5 times more than the control. The ratio drops again to about 7 and 2 for the two scenarios, respectively, and even further in the far future beyond 2500. In pmt_ssp126 and pmt_ssp245 the main melting takes place on the shelf of the Kara Sea (Fig. 5), whereas in pmt_ssp585 the main ice loss is on the shelves of the Laptev and East Siberian seas. The geographic patterns of the melting are essentially constant over time (not shown). In agreement with earlier findings (Overduin et al., 2012, 2019; Angelopoulos et al., 2019), the largest melting is found near the coast, where the water is shallowest and the submerging took place most recently. Thawing in these regions was initiated comparatively recently, and thus the ice-bonded permafrost table (IBPT) is still close to the sediment–water interface.

Figure 4 Melting rate ratios relative to pmt_pre for the main scenarios. The 10-year averages. Grey shading marks the historical period.

Figure 5 SSPF ice column melting for the three main experiments in the period 1850–3000. Green dots mark points entirely without SSPF ice, and blue dots mark points where all SSPF ice melted before the year 3000.

By the year 3000 the initial area of SSPF (2.89×10⁶ km²) is reduced to 2.63×10⁶, 2.55×10⁶, 2.37×10⁶, 2.22×10⁶, and 1.65×10⁶ km² for the two control runs and the three SSPs, respectively. The disproportionally large loss of SSPF area between 1850 and 2100 (Table 2) compared to the volume loss is explained by the rather large area of thin SSPF ice (0.77×10⁶ km² has <10 m of ice). In general, the relative loss of SSPF area is larger than that of SSPF volume and most pronounced for the intermediate scenario experiments pmt_ssp126 and pmt_ssp245. The two control experiments (pmt_pre and pmt_curr) apply a cyclic forcing which means that there are no long-term geographical shifts in the benthic temperatures. This leads to an essentially fixed area to volume melting ratio.

Table 2 Overview over volume (left) and area (right) of SSPF ice in the main and control experiments at selected times. Numbers in brackets indicate relative loss compared to 1850.

Initially the highest SSPF ice saturation (≈86 %, Fig. 1b) is found in a layer at around 6 m depth. In all main experiments the ice saturation of this layer is reduced from about 2050 onward (Fig. 6) but to a very different degree. In pmt_ssp1-2.6 and pmt_ssp2-4.5 it is reduced to a saturation of ≈55 % and ≈30 % in the year 3000, respectively, whereas in pmt_ssp5-8.5 the ice at this depth is completely melted shortly after the year 2100. This layer effectively insulates the deeper ice and thus prevents ice melting in the deeper layers. In pmt_ssp585 the IBPT quickly deepens after the year 2100 and reaches a depth of ≈100 m before the year 2500. The melting rate of SSPF ice however decreases during this period due to the lower concentrations of SSPF ice between 6 and 100 m depth.

Figure 6 Vertical distribution of SSPF ice saturation (main panels) and sea ice concentration (top panels) for the scenario experiments. Yearly averages for the 242 points which contained ice in the initial conditions (1850). “Depth” is measured downward from the sea floor.

SSPF ice is melted from below by geothermal heat in addition to the climatically driven melting from above. All ice in the lowest layer (layer 22, centre at 850 m depth, initially 0.69 % saturated) is gone by the year 3000, and in layer 21 (centre depth 550 m) the ice saturation decreases from an initial 14.5 % to 12.2 % by the year 3000. Effects of geothermal heat flux are also seen in the above layer(s), but in these layers the melting is caused by a combination of heat from below (geothermal) and above (climate).

Figure 7 SSPF ice saturation (a, lat–depth planes), sea ice concentration (a, lat–time planes), sediment temperature change (b, lat–depth planes), and atmospheric 2 m temperature change (b, lat–time planes) along a section across the East Siberian Shelf (153^∘ E, shown on the insert in a) for different time periods. In (b) time slice at 1900 represents the 1850–1900 change, and subsequent slices represent the changes relative to the slice before. Black lines at constant depth are model layer boundaries. Magenta line shows the water depth along the section. Shown are (differences between) yearly averages.

The development of melting and temperature for a selected section across the East Siberian Shelf in pmt_ssp585 is illustrated in detail in Fig. 7, where seven time slices of (panel a) SSPF ice saturation and sea ice concentration (SIC) and (panel b) sediment and air temperature changes are shown. Before the year 2200 salinity-induced taliks (unfrozen lenses) are found in the layers centred at 38 and 150 m. By the year 2100 only a little ice at 2–3 m depth is melted away and the temperature at the same depths rises by ≈2 K preferentially towards the coast. In the 22nd century a dramatic rise in temperature of >8 K in the sediments above ≈5 m is associated with melting of all ice above 10 m. As long as sea ice is present (until about 2100) the temperature rise in the upper sediments is decoupled from the atmospheric temperature rise. After the disappearance of the sea ice, the temperature rise of upper sediments and atmosphere are closely linked together.

The insulating, and thus SSPF preserving, effect of sea ice can also be observed on the pan-Arctic scale (see Fig. S10 and compare Figs. 2, S2, and S3). In the 19th and early 20th centuries the annual mean SIC of the SSPF area is 74 % (68 % of the modelled area, Fig. S2), rapidly dropping in the first half of the 21st century. In SSP1-2.6 this is followed by a slight rise until about 2200 and a stabilization at >50 %. SSP2-4.5 also stabilizes though below 40 % after a drop at a lower rate until about 2150. In SSP5.8-5 the sea ice continues to shrink at a fast rate and is entirely gone at the end of the 21st century except from seasonal ice in the coastal waters of the Laptev Sea. Thus, from this point in time also at pan-Arctic scale, the benthic temperature (Fig. 2) starts to adapt to the higher atmospheric temperatures (Fig. S3).

Figure 8 Temperature difference between ocean bottom (“benthic” temperature) and the uppermost subsea sediment layer as a function of fractional local sea ice concentration for the scenario experiments. Logarithmic density plot (temperature difference in 0.05 K bins, sea ice cover in 0.25 % bins) of monthly mean values. The colours are combined RGB-wise according to the colour triangle. The colour saturation is scaled such that 0 points in a bin correspond to the value 0, and ≥100 points correspond to the value 1. Thus for instance a red point contains an equal number of points from pmt_ssp126 and pmt_ssp245 but none from pmt_ssp585. Black indicates at least 100 points in all three scenarios. Each scenario contributes more than 460 000 time–space points.

The melting of SSPF ice from above is determined by the energy input into the sediments. A quite direct measure for this energy input is the difference between the benthic temperature and the temperature of the upper sediment layer (ΔT). ΔT shows clear patterns when plotted against the local SIC (in time and space; Fig. 8). The local SIC sets a limit for ΔT, which is always essentially 0 K when SIC >70 %. Above this limit, the insulating effect of the sea ice seems to be sufficient to prevent rising temperatures in the water column, and thus the benthic water is close to thermodynamical equilibrium with the upper sediment. With lower SIC the maximum “allowed” ΔT rises approximately linearly, independently of the chosen scenario. The most (ΔT,SIC) points with low SIC are found in pmt_ssp585, but the relation for the maximum ΔT seems to be scenario-independent. This “universal” relation may explain the much larger SSPF melting in pmt_ssp585 compared to the other scenarios.

4 Discussion

4.1 Relation between SSPF ice melting and sea ice

The rate of SSPF ice melting from above is determined largely by the benthic temperature (Fig. 8). The benthic temperature in turn is limited by the presence of sea ice which suppresses energy input into the ocean. Thus a direct relation can be found between the open-sea season and the melting of SSPF ice (Figs. S9, S10, and S11). A similar link was suggested by Dmitrenko et al. (2011) based on the recent observed rise in the benthic temperatures on the East Siberian Shelf. Figure S11 also illustrates the decreasing melting rates from the end of the 21st century due to the complete disappearance of SSPF ice near the benthic surface. In our model the link between the absence of sea ice and the melting of SSPF ice may be exaggerated by the coarse resolution of MPI-OM, implying a zero-curtain effect on the sea-ice–oceanic temperature relation on a relatively large scale. As long as a model grid cell has SIC >0 all energy which goes into the ocean in this grid cell is used to melt the remaining sea ice, preventing a rise in the water temperature. In reality this is very dependent on the heterogeneity of the sea ice. If the remaining sea ice is collected in a corner of a cell, the water may be heated to above freezing point temperatures in the rest of the grid cell. Also the coarse MPI-OM resolution may underestimate the oceanic advection on the Arctic shelf, exaggerating the influence of sea ice on the SSPF on the local scale. Despite these model artefacts, the basic physics of the causal relation chain from disappearance of sea ice to melting SSPF ice seems plausible.

4.2 Model limitations and assumptions

We are ignoring the thermal coupling between sea surface and the sediment at the sea bottom caused by bottom-fast ice (Nicolsky et al., 2012; Osterkamp et al., 1989) since such effects are, due to the coarse resolution of the ESM set-up of JSBACH (Sect. 2.4), never relevant for a significant part of a grid cell.

Freezing of the upper part of subsea sediment below shallow waters is in reality possible (Osterkamp et al., 1989), though this is mainly observed where bottom-fast ice is present. Due to the lack of diffusion of salt in the sediments in our model, we made the assumption that freezing of the below-sea sediments is not happening since in reality the salinity of the porewater of the upper sediments would be in equilibrium with the liquid seawater above. The validity of this assumption decreases with increasing depth into the sediments due to the larger time lag of the salinity entrainment caused by the difference between temperature and salinity diffusivities (Harrison and Osterkamp, 1978). We tested the consequences of this assumption by a series of sensitivity experiments (Table 1), allowing the entire sediment column (pmt_freeze and pmt_freeze126) or the sediment below 30 cm (from layer 4 downward) to freeze (pmt_fr3). These experiments show that freezing occurs mainly in the upper 2–3 layers (not shown). The water volume frozen (years 1850–3000) is 41, 170, and 12 km³ for the three experiments, respectively. This is very small compared to the total volume of SSPF ice (≈ 90–150 × 10³ km³). Due to the energy buffering effect of the freezing–melting cycle in the upper layers, a much larger volume of SSPF ice is preserved: 2.5×10³ km³ less SSPF ice melts by the year 2000 (Fig. S7). For the strong forcing experiments (pmt_freeze and pmt_fr3) the difference to pmt_ssp585 is decreasing after the year 2000 and is almost gone after the year 2300. This “catching up” of the melting is facilitated by the rising benthic temperatures which essentially prevent any freezing of sediments in the later part of these experiments. In the low forcing experiment, pmt_freeze126, freezing takes place at an essentially constant rate throughout the experiments, preserving the energy buffering in the upper layers. This results in a slight but persistently lower rate of SSPF ice melting than in pmt_ssp126. In total, the difference to the main experiments is at most ≈2 % of the total SSPF ice volume. Based on these sensitivity experiments, we conclude that the influence of the “no refreeze” assumption is minor and does not alter our main conclusions. The “all-can-refreeze” and the main experiment bound the realistic melting from below and above, respectively.

In this study, as well as in the study by Overduin et al. (2019), the effects of increasing hydrostatic pressure with increasing depth on the freezing point were ignored. A preliminary estimate suggests that these effects are almost an order of magnitude smaller (0.35 K at 500 m, not shown) than the freezing point depression by salinity. Most of the SSPF ice is found above 500 m, so the effects would be even smaller. With this reasoning Overduin et al. (2019) did not include these effects, and including them here would thus also increase the inter-model inconsistency. Including a freezing point depression by hydrostatic pressure would likely introduce an enhanced initial melting of SSPF ice which would in parts be compensated by lower melting rates later on since sediment temperatures would be lower due to additional energy needed for the additional initial ice melting.

In our set-up, we suppress water advection in the unfrozen parts of the sediments. This is however unlikely to make any difference for our main focus on the thawing from above since the main effect is a heating from above which would tend to stabilize the unfrozen sediment water column in the sediments. An active hydrology may however enhance the influence of the geothermal heating by increasing the effective upward heat flux.

On the Arctic shelf, a broad range of sedimentation rates have been reported. From previous works, Overduin et al. (2019) report rates from near zero to 700 cm kyr⁻¹. For their model, they chose 30 cm kyr⁻¹ for submerged points. Over the experimental time of the present study, this would lead to a sedimentation of ≈34.5 cm of additional sediment or 12 cm during the first 400 years of the experiments when the most interesting dynamics take place. Compared to a typical depth of ice melting of at least 8 m (Fig. 6) and the typical depth to which the seasonal temperature changes penetrate (30 % of surface amplitude at ≈4 m, Fig. S13), the added sediment is thus small, and we therefore neglected sedimentation in the present study. The penetration depth is remarkably constant given the large temperature range of seasonal amplitudes found at different locations (Fig. S12).

The thermal properties of the (dry) sediments were, for consistency with other studies using JSBACH and MPI-ESM, kept equal to the properties of the standard JSBACH “bedrock” (conductivity 2 W km⁻¹, capacity 2 × 10⁶ J km⁻³). These values may not be adequate for subsea sediments, so – inspired by the measurements of Goto et al. (2017) in muddy sediments of Japan – we conducted two additional experiments (pmt_pre_lc and pmt_ssp585_lc), setting the heat conductivity for (dry) sediments to 0.8 W km⁻¹ and their heat capacity to 3.6 × 10⁶ J km⁻³. Adding additional pore water to these values, the sediment thermal properties in these experiments are quite close to those of pure water. Thus these experiments can be regarded as an extreme case for the lowest possible heat entrainment into the sediments. The results (Figs. S8 and S9, equivalent to Figs. 3 and 4, respectively) reveal that the changed sediment properties lead to a much lower melting of SSPF ice. Until the year 2500 only 38 % and 57 % are melted compared to the respective main experiments (pmt_pre and pmt_ssp585, respectively). Thus indeed less SSPF ice will melt and less carbon be released using these parameters. However, the relative ratio between the main and the control experiments (Fig. S9) even increases compared to the experiments using standard sediment thermal properties, emphasizing the influence of climate change. What remains largely unchanged is the time at which the melting starts to accelerate in the high-emission scenario. Also in pmt_ssp585_lc this is around 2100, leading to extreme melting rates in the 22nd century ≈2.5 times those found in the moderate forcing scenario pmt_ssp245 which used the standard sediment thermal properties.

4.3 Area of SSPF

Our present-day (2020) SSPF ice area, $\approx \mathrm{2.74}\times {\mathrm{10}}^{\mathrm{6}}$ km², is larger than the SSPF area reported by Overduin et al. (2019), 2.48×10⁶ km², which in turn is somewhat larger than the ≈ 2×10⁶ km² presented in Sayedi et al. (2020). Part of this difference may be explained by Overduin et al. (2019) presenting a pre-industrial estimate, while Sayedi et al. (2020) refer to the present day. In our model ≈ 0.15×10⁶ km² of SSPF ice area disappears during the 170-year period in between. A total of 76 % of this retreat is also found in the control experiment pmt_pre, indicating that about $\approx \mathrm{0.036}\times {\mathrm{10}}^{\mathrm{6}}$ km² of SSPF ice melts away due to climatic changes within the historical period. As discussed in Overduin et al. (2019) also the idealized benthic temperatures used for the SuPerMAP runs are likely to be too cold in regions of oceanic inflow into the Arctic Ocean such as downstream of the Gulf Stream current system and in the vicinity of the Bering Strait. This may result in larger areas with SSPF in both their and our study than observationally based studies. Both Overduin et al. (2019) and Sayedi et al. (2020) use the more common permafrost definition, based on continuous years of sub-zero temperatures. If we apply this method, our SSPF area is even larger: 3.31×10⁶ km² in 1850 (3.18×10⁶ km² in 2020). That our SSPF area estimates are larger than the previous studies is likely a result of the coarse T63 resolution causing extrapolation of SSPF to locations further offshore (by up to ≈90 km) than included in the model of Overduin et al. (2019). A 90 km wide stripe extending around the pole at 75^∘ N has an area of 0.93 × 10⁶ km². Within this margin, our areas are in agreement with the other estimates.

4.4 Geothermal heat flux and thawing rates (rate of change in IBPT)

The two sensitivity experiments excluding the geothermal heat flux (pmt_pre_0 and pmt_ssp585_0) melt 4.7×10³ km³ and 4.9×10³ km³ (4.0 and 4.2 km³ a⁻¹), respectively, less than the experiments with geothermal heat flux (Fig. S6). For ssp_pre this is more than half of the total melting, indicating that geothermal heat flux is more important than climate forcing for the melting. These two experiments represent the extremes of the forcing, and since the geothermal melting is virtually identical, it can be concluded that climate-caused melting from above and geothermally caused melting from below are largely independent from each other and thus can be treated independently.

The rate of change in the IBPT controls the rate at which carbon become available for degradation. It can be roughly estimated by $\mathrm{\Delta }\mathrm{IBPT}=\left(T_{\mathrm{total}}-T_{\mathrm{noheat}}\right)/\stackrel{\mathrm{‾}}{\mathrm{por}}$, where T_total and T_noheat are the SSPF ice melting rates in the main experiment and the experiment without geothermal heat flux, respectively, and $\stackrel{\mathrm{‾}}{\mathrm{por}}$ is the average porosity over the layers where SSPF is thawed. Assuming $\stackrel{\mathrm{‾}}{\mathrm{por}}=\mathrm{0.35}$ and ignoring the shrinking of SSPF area with time, the average ΔIBPT for points with SSPF ice is 0.8, 0.9, 1.2, 1.7, and 5.3 cm a⁻¹ for pmt_pre, pmt_curr, pmt_ssp126, pmt_ssp245, and pmt_ssp585, respectively, for the 1850–3000 period. The ΔIBPT for the 1850–2009 period is 1.2 cm a⁻¹ and can be compared to present-day measurements. Overduin et al. (2015) measured thawing rates between 0 and 16 cm a⁻¹ off the Muostakh Island and Overduin et al. (2012) 0.5–1.4 cm a⁻¹ off Barrow. Both these study sites were in shallow waters very close to the coast. The model hindcast study of Angelopoulos et al. (2019) presented thawing rates of 1.2–1.5 cm a⁻¹ and thus very similar to the rates found there. In the 22nd century peak of pmt_ssp585 the thaw rate estimate is 11.4 cm a⁻¹.

5 Conclusions

Here, for the first time, a land component of an ESM was used to project the development of subsea permafrost (SSPF) ice until the year 3000 by forcing it with extended CMIP6 scenarios spanning the likely range of climate change. Using an ESM component has the advantage over more specific SSPF models that it can be directly coupled to the carbon cycle to study climatic feedbacks. Though the inter-scenario differences are minor in the present century, in the 22nd century, the strong forcing scenario SSP5-8.5 strongly diverges by melting all SSPF ice above ≈100 m depth, whereas the moderate and low forcing scenarios still leave SSPF ice below ≈8 m. This highlights the need to go beyond 2100, the usual timescale on which climate projections are made, to get the response of slowly acting climate components such as SSPF to a changing climate. The large loss of SSPF ice in SSP5-8.5 is closely linked to the disappearance of local sea ice. Specifically, the length of the open-water season controls the temperature rise of the ocean water on the Arctic shelf. This temperature rise in turn determines the energy input into the ocean bottom sediments and thus the melting of SSPF ice. Therefore our results suggest that the length of the local open-water season may be used as an easy observable proxy for the rate of melting of SSPF ice.