Permafrost and Gas Hydrate Stability Zone of the Glacial Part of the East-Siberian Shelf

By using thermal mathematical modeling for the time range of 200,000 years ago, the authors have been studying the role the glaciation, covered the De Long Islands and partly the Anjou Islands at the end of Middle Neopleistocene, played in the formation of permafrost and gas hydrates stability zone. For the modeling purpose, we used actual geological borehole cross-sections from the New Siberia Island. The modeling was conducted at geothermal flux densities of 50, 60, and 75 mW/m2 for glacial and extraglacial conditions. Based on the modeling results, the glaciated area is characterized by permafrost thickness of 150–200 m lower than under extraglacial conditions. The lower boundary of the gas hydrate stability zone in the glacial area at 50–60 mW/m2 is located 300 m higher than the same under extraglacial conditions. At 75 mW/m2 in the area of 20–40 m isobaths, open taliks are formed, and the gas hydrate stability zone was destroyed in the middle of the Holocene. The specified conditions and events were being formed in the course of the historical development of the glacial area with a predominance of the marine conditions peculiar to it from the middle of the Middle Neopleistocene.


Introduction
The information related to the submarine cryolithozone of the East Siberian Arctic shelf, namely, its distribution, the table depth, and the frozen deposit thickness, is essential to drawing up the forecasting scenarios of climate warming. Degrading submarine frozen deposits are considered as the greenhouse gas source of a planet-importance level. Without knowing cryolithozone expansion and structure, it is impossible to assess the prospects for efficient use of the Northern Sea Route, an international transport and logistics system of the nearest future, the environmental footprint of oil and gas resources exploration, and development within the East Arctic shelf.
The drilling methods were applied to study the cryolithozone of the Laptev and East Siberian seas shelf only within the near-shore zone. The maximum distance from boreholes to the coast is 25 km. The drilling depth, as a rule, does not exceed 40-50 m. The boreholes with a depth of 216 and 203 m (at the sea depths of 12 and 14 m, respectively) did not penetrate the permafrost bottom [1]. Thus far, no boreholes have been drilled to penetrate the entire thickness of permafrost. The seismic method's application did not leave the stage of method development. The electrical exploration allowed for revealing a vast talik area in the Khatanga Bay and Nordvik Bay [2]. The lack of the possibility to obtain zone under the influence of glaciation and subsequent sea transgression. According to the results presented, the upper boundary of the methane hydrates stability zone in Yamal could have reached the ground surface during the periods of glacial conditions. This brief review shows that some models presenting the current state of the cryolithozone of the East Siberian sector within the Arctic seas have been published so far. With all their advantages, they have certain downsides. The downsides are determined by insufficient coverage of studies of the submarine cryolithozone, shelf geology and paleogeography. The first downside is relevant for all studies currently existing. We mean hypothetical, schematic, and often groundless character of the given geological sections serving as the basis for the calculations. The second downside is the lack of testing for the paleogeographic scenario and specified deposit properties. However, testing serves as justification for further modeling performed on its basis. These circumstances cause doubts about the obtained results in terms of their feasibility. It is often impossible to verify them due to unavailability of relevant actual data. Finally, the third downside: all the models are based on the idea of no Middle Neopleistocene glaciation in the East Siberian sector of the Arctic. However, the information on the glacial genesis of ground ice on the New Siberia island was published in 2003 [35], considered as tentative, and in 2006-2008 [36,37] as established. B.G. Golionko et al. [38] proved the glacial nature of deformations in the Cretaceous-Middle Neopleistocene strata of the New Siberia and Faddeevsky Islands.
Based on the foregoing, the purpose of the research is stated as follows: the assessment of evolution and current distribution of permafrost and GHSZ, their upper and bottom boundaries in the glacial part of the East Siberian shelf. The specified assessment is performed using thermal mathematical modeling and provides for comparison of its results with the permafrost and GHSZ parameters under nonglacial conditions. To achieve the determined goal in the glacial area and under nonglacial conditions of the East Siberian shelf, the following tasks were to be solved: (1) To create a scenario of geological development of the area of study. (2) To test the created scenario.

Study Area
The area of study ( Figure 1) is located within the East Siberian shelf. The landscapes of New Siberia Island and Faddeevsky Peninsula are represented by the arctic tundra and the arctic desert. The area is classified as nonglacial (periglacial) because no giant ice sheets were formed here in the Middle-Late Neopleistocene. However, as indicated above, a single local ice sheet occurred. Its center was presumably located on the Zhokhov, Vilkitsky, and Bennetta Islands ( Figure 2) [37]. Modern absolute heights of the earth's surface in the glaciation area vary from +40 to −40 m. The glacier covered the area of 150,000-200,000 km 2 [39] and caused the development of alternating glacial isostatic movements. They resulted in sea transgressions and regressions, which formed the stratum of coastal-marine and marine deposits of the Kanarchakskaya suite. It consists of upper and lower subsuites separated by a horizon with ice ground and sheet ice [40]. The sedimentation depth of the lower subsuite varied from a few meters to 50 m. The roof of the Kanarchak sediments in the north of the Faddeevsky and New Siberia Islands is relief-forming. Its absolute heights are 20-40 m. Based on the biostratigraphic data and the results of 230Th/234U-dating, the glaciation age was determined as the end of the Middle Neopleistocene (MIS-6) [41]. A glacier deformed the Upper Cretaceous-middle Neopleistocene strata of the New Siberia Island and Faddeevsky Peninsula [38]. The sediments of the upper Kanarchakskaya subsuite overlap the deformed stratum.

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: Geosciences 2020, 10, x FOR PEER REVIEW 5 of 22

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: At the boundary between frozen and thawed deposits, the condition of equality of the ground temperature to the water freezing temperature TF and Stefan condition for the moving boundary of phase transitions are allowed: At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, (3) At the lower boundary of the area calculated (HS = 1500 m), the geothermal flux G is set: at z = HS: λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux.
Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost. Figure 3. Pressure-temperature equilibrium relationship in the phase diagram of the water-icemethane-hydrate system [42].

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: (1) At the boundary between frozen and thawed deposits, the condition of equality of the ground temperature to the water freezing temperature TF and Stefan condition for the moving boundary of phase transitions are allowed: At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, (3) At the lower boundary of the area calculated (HS = 1500 m), the geothermal flux G is set: at z = HS: λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost. Figure 3. Pressure-temperature equilibrium relationship in the phase diagram of the water-icemethane-hydrate system [42].

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: (1) At the boundary between frozen and thawed deposits, the condition of equality of the ground temperature to the water freezing temperature TF and Stefan condition for the moving boundary of phase transitions are allowed: At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, (3) At the lower boundary of the area calculated (HS = 1500 m), the geothermal flux G is set: at z = HS: λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost. Figure 3. Pressure-temperature equilibrium relationship in the phase diagram of the water-icemethane-hydrate system [42].

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: Сi At the boundary between frozen and thawed deposits, the condition of equality of the ground temperature to the water freezing temperature TF and Stefan condition for the moving boundary of phase transitions are allowed: At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, (3) At the lower boundary of the area calculated (HS = 1500 m), the geothermal flux G is set: at z = HS: λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost. Figure 3. Pressure-temperature equilibrium relationship in the phase diagram of the water-icemethane-hydrate system [42].

Methods
The research was conducted using numerical mathematical modeling. It also provides f creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we u a one-dimensional model of thermal processes in bottom sediments taking into account the pha transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional he equation: Сi × ∂ T/∂t = (∂/∂z)(λi × ∂T/∂z), ( At the boundary between frozen and thawed deposits, the condition of equality of the groun temperature to the water freezing temperature TF and Stefan condition for the moving boundary phase transitions are allowed: λth × ∂T/∂z = G, ( The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сi volumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes o of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moistu L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-d bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the dep and composition of sediments and are set in accordance with the accepted geological model (s Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of t thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state the GHSZ within the permafrost area strongly depends on the ground temperature and pressu that is why the calculations were conducted with the phase curves of methane hydrate under t corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer w calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as t lithostatic pressure at the bottom of permafrost. Figure 3. Pressure-temperature equilibrium relationship in the phase diagram of the water-icemethane-hydrate system [42].
At the boundary between frozen and thawed deposits, the condition of equality of the ground temperature to the water freezing temperature T F and Stefan condition for the moving boundary of phase transitions are allowed: at z = Z F : λ th × ( Geosciences 2020, 10, x FOR PEER REVIEW 5 of 22

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: Сi At the boundary between frozen and thawed deposits, the condition of equality of the ground temperature to the water freezing temperature TF and Stefan condition for the moving boundary of phase transitions are allowed: At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, (3) At the lower boundary of the area calculated (HS = 1500 m), the geothermal flux G is set: at z = HS: λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost.

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: Сi At the boundary between frozen and thawed deposits, the condition of equality of the ground temperature to the water freezing temperature TF and Stefan condition for the moving boundary of phase transitions are allowed: At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, (3) At the lower boundary of the area calculated (HS = 1500 m), the geothermal flux G is set: at z = HS: λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost.

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: Сi At the boundary between frozen and thawed deposits, the condition of equality of the ground temperature to the water freezing temperature TF and Stefan condition for the moving boundary of phase transitions are allowed: At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, (3) At the lower boundary of the area calculated (HS = 1500 m), the geothermal flux G is set: at z = HS: λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost.

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: Сi At the boundary between frozen and thawed deposits, the condition of equality of the ground temperature to the water freezing temperature TF and Stefan condition for the moving boundary of phase transitions are allowed: At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, (3) At the lower boundary of the area calculated (HS = 1500 m), the geothermal flux G is set: at z = HS: λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost.

Methods
The research was conducted using numerical mathematical modeling. It creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for ga a one-dimensional model of thermal processes in bottom sediments taking into transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one equation: Сi × ∂ T/∂t = (∂/∂z)(λi × ∂T/∂z), At the boundary between frozen and thawed deposits, the condition of equa temperature to the water freezing temperature TF and Stefan condition for the m phase transitions are allowed: λth × ∂T/∂z = G, The notations used: T-ground temperature; t-time; z-depth from the bo volumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; su of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volum L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × bulk density; TB-temperature on the bottom sediments surface; G-the geotherm Physical and thermal properties of sediments used for numerical modeling de and composition of sediments and are set in accordance with the accepted geo Section 3.2). We assume that the sediments pores are completely filled with water The calculations of the thermal state of bottom sediments are accompanied thermal dynamic boundaries of the methane gas hydrates stability zone (Figure 3 the GHSZ within the permafrost area strongly depends on the ground tempera that is why the calculations were conducted with the phase curves of methane corresponding lithostatic pressure. In accordance with [43], the pressure in the calculated as the hydrostatic column pressure beneath the frozen layer at the stud lithostatic pressure at the bottom of permafrost.

Methods
The research was conducted using numerical mathematical modelin creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions a one-dimensional model of thermal processes in bottom sediments taking transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the equation: Сi × ∂ T/∂t = (∂/∂z)(λi × ∂T/∂z), At the boundary between frozen and thawed deposits, the condition of temperature to the water freezing temperature TF and Stefan condition for t phase transitions are allowed: At the upper boundary of bottom sediments corresponding to the botto at z = 0: T = TB, At the lower boundary of the area calculated (HS = 1500 m), the geotherm at z = HS: λth × ∂T/∂z = G, The notations used: T-ground temperature; t-time; z-depth from t volumetric heat capacity of the ground; λi-soil thermal conductivity coefficie of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative vo L-specific heat of freezing and melting of water in the ground pores (L = 3 bulk density; TB-temperature on the bottom sediments surface; G-the geot Physical and thermal properties of sediments used for numerical modeli and composition of sediments and are set in accordance with the accepted Section 3.2). We assume that the sediments pores are completely filled with w The calculations of the thermal state of bottom sediments are accompan thermal dynamic boundaries of the methane gas hydrates stability zone (Fig the GHSZ within the permafrost area strongly depends on the ground tem that is why the calculations were conducted with the phase curves of meth corresponding lithostatic pressure. In accordance with [43], the pressure in calculated as the hydrostatic column pressure beneath the frozen layer at the lithostatic pressure at the bottom of permafrost.
At the upper boundary of bottom sediments corresponding to the bottom surface: At the lower boundary of the area calculated (H S = 1500 m), the geothermal flux G is set: at z = H S : λ th × Geosciences 2020, 10, x FOR PEER REVIEW 5 of 22

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: Сi At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone (Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost.

Methods
The research was conducted using numerical mathematical modeling. It also provides for creation of a scenario of geological development of the region.

Description of the Numerical Model
To calculate the thermal state of permafrost and thermobaric conditions for gas hydrates, we use a one-dimensional model of thermal processes in bottom sediments taking into account the phase transitions between frozen and thawed ground [18,19].
Heat distribution in bottom sediments is described by solving the one-dimensional heat equation: Сi At the upper boundary of bottom sediments corresponding to the bottom surface: at z = 0: T = TB, λth × ∂T/∂z = G, (4) The notations used: T-ground temperature; t-time; z-depth from the bottom surface; Сivolumetric heat capacity of the ground; λi-soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); WTot-relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρD-dry bulk density; TB-temperature on the bottom sediments surface; G-the geothermal flux. Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone (Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost.
The notations used: T-ground temperature; t-time; z-depth from the bottom surface; C i -volumetric heat capacity of the ground; λ i -soil thermal conductivity coefficient; subscript i takes one of the values "f" (frozen ground) or "th" (thawed ground); W Tot -relative volumetric soil moisture; L-specific heat of freezing and melting of water in the ground pores (L = 3.35 × 10 5 J/kg); ρ D -dry bulk density; T B -temperature on the bottom sediments surface; G-the geothermal flux.
Physical and thermal properties of sediments used for numerical modeling depend on the depth and composition of sediments and are set in accordance with the accepted geological model (see Section 3.2). We assume that the sediments pores are completely filled with water or ice.
The calculations of the thermal state of bottom sediments are accompanied by estimates of the thermal dynamic boundaries of the methane gas hydrates stability zone ( Figure 3) [42]. The state of the GHSZ within the permafrost area strongly depends on the ground temperature and pressure, that is why the calculations were conducted with the phase curves of methane hydrate under the corresponding lithostatic pressure. In accordance with [43], the pressure in the ground layer was calculated as the hydrostatic column pressure beneath the frozen layer at the study depth and as the lithostatic pressure at the bottom of permafrost.
The boundary condition on the surface of the bottom sediments (3) is determined by the transgressions-regressions periods with due account for the sea level changes, as well as by glacial conditions having existed over the last 200,000 years. Curves of temperature changes on the surface of the bottom sediments, used for modeling submarine permafrost, are shown below in Section 4.2. The geothermal flux intensity was taken equal to G = 50, 60, or 75 mW/m 2 , depending on the numerical calculation, which corresponds to the values of the heat flux for this region [44,45].
The numerical model implementation (1)-(4) is based on the sweep method on a discrete calculation grid with a vertical step of 0.5 m and an implicit time scheme with a step of 1 month.

Methods for Generating a Scenario of Geological Development for Research Area
The scenario of the geological development consists of a scenario of sea level fluctuations, paleogeographic (paleotemperature) scenario, and a geological model for the research area.
To set the geological structure and physical and thermophysical properties of deposits in glacial and extraglacial conditions, the data on the geological structure of the New Siberia Island were used.
The sea level scenario is built based on glacioeustatic curves of the World Ocean level fluctuations and their transformation following the regional data [11]. This approach has been applied in mathematical modeling of the evolution and the current state of the Laptev and East Siberian seas cryolithozone since 2003 [7,9]. In the present studies, we used the World Ocean level fluctuations curve and the regional data.
The eastern sector of the Russian Arctic has no paleotemperature data for several time intervals. This circumstance precludes from building a scenario of air or ground temperature dynamics using only regional data. For this purpose, to build such a scenario, the method with the isotopic paleotemperature curves obtained from the East Antarctic and Greenland glacial cores reflecting the course of global climate fluctuations, is used. Such curves are to be transformed into a regional paleotemperature scenario (model) with the help of regional paleotemperature reconstructions used as paleotemperature benchmarks. The isotopic paleotemperature curves allow for building continuous in time regional curves of the air and ground temperature dynamics using the discrete regional data [10,11]. As the latter, we used the paleotemperature reconstructions for the coastal lowlands of the north of Yakutia and the New Siberian Islands. They are created according to the width of elementary ice veins in the ice wedges [46], according to the ratio of the main rock-forming minerals [47], according to the isotopic composition of ice wedges [48][49][50], according to geothermal observations in deep boreholes [51], according to a set of paleoecological data [52,53].
Such a scenario includes a family of paleotemperature curves. We construct paletemperature curves for different isobaths and regions based on the shelf zoning in terms of bathymetry, tectonic structures, and geological development [10,11,54]. This method based on transformation of the isotope curve of Vostok station (East Antarctic), was widely used earlier [6,9,16] and is being used in a modified form now [19,23,25].

Geological Model
Assessment of the impact glaciation has on the current parameters of the permafrost and GHSZ and their evolution makes it necessary conduct such assessment under otherwise identical conditions. Therefore, to set the geological structure, sediments stratification, composition, physical and thermophysical properties of deposits in glacial and extraglacial conditions, the data on the geological structure of New Siberia Island were used. The density of the geothermal flux following the map [44], as indicated above, was 50, 60, and 75 mW/m 2 .
The area under research falls within the young Epikimmerian platform. Its pre-Aptian basement within this territory includes the Middle Carboniferous-Valanginian structural stage (C 2 -K 1 v). The platform cover is of Aptian-Cenozoic age (K 1 ap-KZ) and has the maximum thickness in depressions. Within uplifts, some horizons decrease in thickness or disappear from section [55]. New Siberian depression includes the New Siberia Island, almost the entire Faddeevsky Island, as well as the entire water area to the east of the Lyakhovsky Islands. New Siberian depression comprises terrigenous strata of the Aptian-Albian and Late Cretaceous-Cenozoic age with a thickness of up to 5 km and more. At the base of the De Long uplift, which includes the De Long Archipelago Islands, a site of the Epicarelian plate and a folded system of the Mesozoic time exist [55].
The plate cover (Late Cretaceous-Cenozoic structural floor) includes deposits of the Late Cretaceous (Bunginskaya and Derevyanogorskaya suites), Paleogene-Miocene (Anzhuyskaya suite), as well as polygenetic sediments of the Pliocene-Quaternary age [56]. The uppermost part of the section is composed of marine sediments of the Nerpichinskaya (Lower-Middle Neopleistocene) and Kanarchakskaya suites (Middle-Late Neopleistocene) [40,41]. When calculated, to ensure the lower boundary of the calculated area does not influence the permafrost bottom dynamics, it was assumed that the maximum vertically set area of the geological model should be 1.5 km. This means the properties of the deposits are to be set down to a depth of 1.5 km. According to V.K. Dorofeev et al. [56], the plate cover thickness varies from a few meters to 300 m onshore and in coastal zones to a few kilometers outside the shelf (sea depths over 500 m). The calculations take the data for two boreholes, which were subject to temperature logging and the estimation of permafrost thickness. This is borehole c-A with a depth of 200 m. It is located on the southern coast of New Siberia Island in the Derevyannye Gory tract (Figures 1 and 4). Based on this borehole the Upper Cretaceous sediments thickness is more than 200 m. The deposits of the Derevyanogorskaya suite (K 2 dr) occurring from the surface with a thickness of 130 m are represented by tuffaceous sands, silts, brown coals with interbedded clay. The total thickness of brown coals is about 40 m. Below this suite, the Bunginskaya suite rocks (K 2 bn), represented by clays, mudstones, silts, and sands, were exposed by drilling. It is assumed that the section below is represented by sandstones, mudstones, silts of the Jurassic age (by analogy with the neighboring areas).
In the area of Cape of Rozhin (the extreme southwestern tip of the New Siberia Island), borehole c-B (borehole depth is 80 m) penetrated modern sea sands and silts with 5 m thickness. The Kanarchakskaya suite silts with a total thickness of 70 m, at 59 m depth containing glacial ice of 40 cm thickness underlie them. Below the sands of the Nerpichinskaya suite were penetrated. Their thickness is set equal to 60 m. Below gravels and sands of the Anzhuyskaya suite are located. Deposits of the Derevnogorskaya and Bunginskaya suites underlie them. Jurassic sandstones, mudstones, and silts are represented from a depth of 500 m and more ( Figure 4). The lower part of the given section is compiled based on the materials of the regional generalization [55]. The solid line of the temperature curve is plotted according to the results of geothermal observations in the borehole, the dashed line-following the geothermal gradient in the depth interval of 50-80 m.
Thermophysical characteristics (Table A1) were assigned based on the lithological composition, geological and genetic type and age of the sediments according to the data obtained within the adjacent territories [57][58][59][60] regarding the known physical characteristics of the deposits from boreholes c-A and c-B.

Sea Level Dynamics Scenario
The scenario of sea level dynamics in extraglacial conditions. For this research ( Figure 5), we took the curve of the World Ocean level fluctuations as the initial one [61]. We faced an extreme shortage of regional data and their inequality to transform this curve into the curve of the East Siberian Sea. Only the course of the Flandrian transgression in the Laptev Sea has been sufficiently completely characterized by such data [62]. Particularly, we mean the information about the dating and altitude of the points of continental sedimentation change into marine during the Flandrian transgression. To determine the shelf drainage limit within cryochron MIS-2, we used the materials of permafrost interpretation of the seismoacoustic data on the shelf edge [63].
The information on the level fluctuations within MIS-3 is contradictory. Only faunistic data in the East Siberian sector of the Arctic are considered as reliable [64]. They show that the level did not exceed abs. heights of −40 m (4 in Figure 5) Contrast alternating level fluctuations during MIS-5 are provided with only one date (IR-OSL 111 ± 5 thous. years ago) [65]. It corresponds to its low state for MIS-5d within the western part of the Laptev Sea (4 in Figure 5). The history of the sea level fluctuations at the end of the Middle Neopleistocene has been also poorly characterized by the regional data. When characterizing the sea level in the last interglacial of the Middle Neopleistocene, it was taken into account that the data on the Konkovskaya transgression in the lower reaches of the river Kolyma are insufficiently substantiated in terms of age [66]. The use of the listed data made it possible to build a scenario variant (curve) for the extraglacial conditions (1 in Figure 5). The scenario of the earth surface and sea level dynamics in the glacial area. Glacioisostatic earth crust movements accompany glaciation. During glaciation, the glacier bed sinks under glacial load. When the glacier degrades and the glacial load decreases, the bed rises. The processes in the upper mantle (asthenosphere) cause these alternating movements. It has the ability of plastic flow within a gravity field [67]. Within the glacial area, the glacioisostatic curve is considered as one of the main parameters to reconstruct the hypsometry of the glacial bed, shores, and their position relative to the sea level. It characterizes the course of glacioisostatic movements in time, that is, the progress of changes in the glacial bed and shores hypsometry [54].
For the glacioisostatic curve plotting, the main thing is reconstruction of the altitude position of the glacial bed in the pessium MIS-6. As the initial data for reconstruction we took correlation of ice sheet with the deepest part of the section of the lower subsuite of the Kanarchakskaya suite [40], as well as the sea depth (50 m) according to the malacofauna-related data [55] and the shelf drainage limits in MIS-6. As far as the East Siberian shelf is concerned, the latter were taken following the data on shelf drainage in MIS-2 [62] and already mentioned results of the permafrost interpretation of the seismoacoustic data on the edge of the Laptev Sea shelf (120 m) [5]. Based on the accepted data the abs. height of the glacial bed in MIS-6 was minus 170 m. Figure 5 shows two glacioisostatic curves. The following data were used to plot them. Particularly, the abs. heights of the glacial bed in MIS-6 and the bottom of sediments comprising ice sheet, as well as the 14 C date of mammal bone residues of the mammoth complex on the New Siberia and Faddeevsky Islands. Based on 14 C data, marine sedimentation not later than 54 ka was replaced with accumulation of continental sediments of the Late Neopleistocene Ice Complex [40]. In terms of calendar time, this date corresponds to the age of at least 59,000-60,000 years ago [68]. At this time, the shelf was drained to 80-100 m isobaths ( Figure 5). The beginning of the regression accrues to the time of 75-70 ka. Therefore, the islands in their current form could have been formed from 75-70 to 60 ka BP. The abs. heights of the roof of the lower subsuite of the Kanarchakskaya suite at the Cape Kamenny are 36-40 m [40]. It underlies sheet ice and was a glacier bed within MIS-6. In other places in the northern part of the islands, this stratum is located lower. Additionally, in the south the New Siberia Island goes below the sea level. The above data and modern abs. heights of the bottom of the lower Kanarchakskaya suite underlying the ice sheets provided the basis for plotting the glacioisostatic curves (2 and 3, Figure 5, for the south and north of the island, respectively). Curve 3 was selected as the main one.
According to curve 3 in Figure 5 and its relationship with the curve of the sea level fluctuations (1 in Figure 5) the glacial area comprised a sea basin in the interval of 180-70 ka. Its level correlated with the glacioisostatic curve. Under extraglacial conditions, subaerial environments were found at this time. Within the interval from 70,000 years ago to the present time according to the ratio of the indicated curves (3 and 1 in Figure 5), the glacial area, as well as the extraglacial area, were subject to the subaerial conditions. The feature of the glacial area was formation (due to glacioisostatic uplift) of a marine terrace up to 40 m high [40].
Thus, the main changes in the course of the sea level fluctuations are described by the ratio of the glacioisostatic curve and the sea level fluctuation curve under extraglacial conditions.

Paleotemperature Scenario
Setting the initial temperature conditions. The part described above shows that during the calculation period the sea conditions existed for a longer time within the glacial area rather than outside this area. The marine periods are the most important warming factor for permafrost formation. To set the initial conditions in the glacial area and outside it, the following should be kept in mind. Within the period preceding the calculated one (200,000 years ago-the present time), the continental regime prevailed in extraglacial conditions [55]. Predominance of the continental regime with cryochron temperature at the level of −23 −27 • C ( Figure 6) determines the long-term existence of submarine permafrost under the extraglacial conditions. According to the data [9], the period of continuous existence of permafrost on the shelf within 0-50 m isobaths was estimated for the last 400,000 years.
The glacial area developed differently. During accumulation of the Nerpichinskaya suite (the first half of the Middle Neopleistocene according to [40]) there was the sea. The sea having existed for hundreds of thousands of years determines lack of permafrost during this long period.
The abovesaid made it possible to set the initial values of the rock temperature equal to −1.8 within the glacial area and −12 • C outside it. The first value roughly corresponds to the current average annual temperature of the bottom water to north of the New Siberia Island, the second value corresponds to the current average annual ground temperatures of the New Siberia Island according to the geothermal observations in boreholes c-A and c-B. The set value for the extraglacial conditions is the highest of possible. It corresponds to the thermochron values. Its duration could hardly exceed 10,000 years. Thus, the initial conditions were set taking into account geological development, both in the glaciated area and outside its boundaries.
Construction of a paleotemperature scenario for the extraglacial area. The scenario allowing for setting the upper boundary conditions when modeling (Figure 6a,c,e) was established based on the methodology in [10,11]. We used the studies, which made it possible to reconstruct the temperature conditions of sedimentation on the New Siberian Islands and the coastal lowlands of northern Yakutia within the different time intervals [46,47,49,50,52,53,65].
The main initial paleodata for MIS-2 were the reconstructed sediment temperatures: −27 • C in the modern Arctic tundra at Cape Svyatoi Nos; −22 −26 • C in woodlands of different parts of the lowlands [47]; −22 −27 in woodlands and subarctic tundra of the lower reaches of the Kolyma [46]. The values of temperatures for warming MIS-5e −2 −3 • C [47] are restored for Cape Svyatoi Nos and the range from 0 to −5 • C for conditions from woodlands to arctic tundra [46]. These authors have given the same values for the early Holocene warming. Changes in air temperatures for the conditions of the southern and northern coasts of the Dm. Laptev Strait over the past 200,000 years are given in [49,50,53,65]. Construction of a paleotemperature scenario for the glacial area. The peculiarities that define the construction of the paleotemperature scenario for the glacial area are the determination of the glacier thickness and setting the ground temperature below the glacier. Glaciation, as it is known [69,70], causes glacial bed subsidence, as it is accompanied by displacement of a subcrustal substance under the glacier load. The density of the subcrustal substance is close to 3.3-3.5 g/cm 3 , and of ice −0.9 g/cm 3 [66]. Therefore, to restore isostatic equilibrium during glaciation, the earth's surface lowered by 2/7 of the average ice thickness. Under deglaciation, the opposite process occurs: the earth's surface rises by the same value. Therefore, having determined the vertical movement value for the glacier bed for the period MIS-6 till present, it is possible to determine the glacier thickness.
The value of this movement is 170 and 200 m (2 and 3 in Figure 5). These values were adopted to estimate the altitude ratio between sea level and coastline. The concern arose from the question if it was possible to use them to estimate the magnitude of the postglacial uplift. In the first approximation, the authors consider it is possible. At the same time, we can also take into account the role of glacial exaration of the apex surface, composed of deposits of Kanarchakskaya suite, and incompleteness of deglaciation. The first factor mentioned can be inferred based on the sea depth in the glacial area. These are the first tens of meters. The estimation of the deglaciation incompleteness effect can add an additional 10 m. This value will result from the melting of the largest ice mass (2/7 of 35 m) on the New Siberia Island. Taking into account the abovementioned terms allows estimating the value of MIS-6 glacier bed elevation as 200-220 m, and the glacier thickness as 700-770 m. These are the maximum values. In the south of the New Siberia Island, the glacier thickness may reach 600-630 m and less. In our calculations, we assumed the glacier thickness equal to 600 m.
To estimate the below glacier ground temperature, we used the value of the geothermal gradient obtained from the geothermal observations in the borehole in the village Tiksi (2 • C/100 m [11,71]). Figure 6b,d,f shows that within cryochron MIS-6 the surface temperature under the glacier was 12 • C higher than within the extraglacial area (−14.5 instead of −26.5 • C).
The glacier is not the only warming factor in the glacial area compared to the extraglacial one. Over the period from 116,000 to 70,000 years ago, the glacial area comprised the sea, while the extraglacial area was subject to the continental conditions with the natural environment of the drained shelf. The ground temperature 20 • C lower than within the glacial area caused active shelf freezing.

Testing of Geological Development Scenario
Mathematical modeling at the sites provided with reliable geocryological data usually performs testing. For the glaciation area, we took the boreholes sites on the New Siberia Island (boreholes c-A and c-B) (Figure 1), provided with the geothermal data. Drilling and thermal logging were conducted during the state geological survey at a scale of 1:200,000 in the 1970s. The geological section, required for solving Stephan's task, in the upper 80 (borehole c-B) and 200 m (borehole c-A) was set by the drilling results, and below-following the geological survey data summarized in [55,56,72]. According to the geothermal observations in borehole c-B (located at the southwestern tip of the New Siberia Island), the average annual temperature was −11.5 • C, the temperature at the bottomhole was −7.6 • C. The permafrost thickness, calculated from the geothermal gradient is 250-255 m [71,73]. The paleotemperature curve and the results of test modeling are shown in Figure 7. Correction of rock properties allowed us to obtain correspondence between the field and estimated data.
We also conducted test modeling using the data from borehole c-B. The average annual ground temperature in the 1970s here was −13.2 • C, the bottomhole temperature −2 • C, the permafrost thickness according to thermal logging and geothermal gradient 235-240 m [71]. Test modeling also showed a coincidence of the calculated data with the field data. The calculations encouraged the conclusion that the generated paleotemperature scenario and geological model are realistic and can be used for mathematical modeling of shelf cryolithozone evolution and current state.

Permafrost Distribution and Thickness
The modeling results for the extraglacial and glacial conditions by the scenario (Figure 6) under the heat flux of 50, 60, and 75 W/m 2 are shown in Figure 8 for areas of 5, 20, and 40 m isobaths. Figure 8 shows the dependence of the depth of the permafrost bottom surface (thickness) on the geothermal flux, sea depth (isobaths), and the geological development history (glacial, extraglacial conditions). The geothermal flux density determines the rate of rocks freezing during shelf drainage and, that is especially important, the rate of permafrost degradation during shelf flooding. The dependence of the permafrost thickness on the sea depth is the dependence on the duration of the period of their aggradation during drainage and degradation during shelf flooding. This period within the glacial area is much shorter than outside this area (Figures 9 and 10). According to the mentioned above, the highest current permafrost thickness should be in the area of 5 m isobaths within extraglacial conditions, and the lowest one in the area of 40 m isobaths of the glacial area under increased heat flux. Indeed, within the first indicated area it is 570 m, within the second one, the relict permafrost has completely degraded.
For this reason, the dependence of the permafrost thickness on the value of the geothermal flux and the depth of the sea (isobaths) is reasonable to compare using the example of the extraglacial conditions, where this dependence is clearly defined. Under the extraglacial conditions, per 5, 20, 40 m isobaths, the permafrost thickness varies from 570 to 500 m at the heat flux of 50 mW/m 2 , from 430 to 370 m at 60 mW/m 2 , and from 240 to 220 m at 75 mW/m 2 .
The dependence of the permafrost thickness on the depth under the extraglacial conditions is much weaker. It gradually decreases as the sea depth increases from 5 to 20 and further to 40 m at each heat flux value of 50, 60, and 75 mW/m 2 . The reduction is expressed in values from 10 to 20-25 m.
The dependence of the permafrost thickness on its location within the glacial area or outside is expressed almost as clearly as for the heat flux. The maximum permafrost thickness is intrinsic to 5 m isobath area-about 340 m under the heat flux of 50 mW/m 2 , 220-230 m less than within extraglacial conditions. With this value of the heat flux at 20 and 40 m isobaths, the permafrost thickness is reduced by 200 m compared to the extraglacial conditions. At 60 mW/m 2 , such reduction amounts to 190 m for 5 m isobath and 170 m for 40 m isobath. At 75 mW/m 2 , relict frozen deposits are found only at 5 m isobath. Their thickness here is almost 100, which is 130,140 m less than within extraglacial conditions. However, on 20 and 40 m isobaths, permafrost thaws. Such a significant transformation of the permafrost strata within the glacial area at 20 and 40 m isobaths at 75 mW/m 2 is stipulated by the following.  When immersed under the sea the frozen deposits are presented as ice bonded permafrost. As the degradation proceeds, along with an increase in the ground temperature, it is transformed into a stratum of plastic-frozen deposits (ice bearing permafrost), where ice is contained in the form of inclusions. Figure 8 shows the results of such degradation, namely, the gradual decrease in the content of ice bonded permafrost in the strata due to an increase in unfrozen water volume and reduction in the ice amount, at the heat flux of 50 and 60 mW/m 2 , both for the glacial and extraglacial conditions. The results of the same degradation, together with temperature increase, change in the unfrozen water and ice ratio, are reflected at the heat flux of 75 mW/m 2 for 5 m isobath within the glacial conditions. The situation is different at 20 and 40 m isobaths within the glacial area. The degradation having sharply accelerated proves that temperature in the strata under degradation increases to the values corresponding to the freezing-thawing point of the pore water, that the phase transitions in it are minimized and that rate of permafrost thawing sharply grows. Rate growing leads to complete permafrost thawing at 20-40 m isobaths.
Thus, with the heat flux of 75 mW/m 2 , discontinuous or sporadic relict permafrost can occur in the glacial area, while within extraglacial conditions, continuous relict permafrost of 200-240 m thickness can occur.

Evolution of Permafrost
The modeling results show that the most significant impact the glaciation had on the evolution of the permafrost (Figures 9 and 10). Under extraglacial conditions, permafrost continuously existed throughout the entire estimated period (200 kyr before present). Within the glacial area, on the contrary, there was a break. Permafrost's existence fell within the intervals of 190-125 and from 70 to 60 ka. During the first interval they formed and existed under the glacier, the second was the interval of their subaerial existence. Its beginning is connected to the sea regression in MIS-4. The age of modern permafrost is reckoned by 60,000-70,000 years.
Permafrost formation and degradation have resulted from the climate and sea level fluctuations. Therefore, the fluctuations are evidenced in the permafrost bottom depth and thickness (Figures 9  and 10). In this case, they are evidenced in a smoothed form and with a delay. Smoothness is most clearly noted at 5 m isobath (Figure 9). The lag of the MIS-5e optimum (125,000 years ago) or the MIS-2 pessium (20,000 years ago) (Figure 6), for example, both in the extraglacial and glacial conditions, is approximately 8000-9000 years. This is in good agreement with the results of other studies [9,74]. The highest permafrost thickness, as evidence of the pessiums MIS-6 and MIS-2, is 450 and 380 m, respectively, within the extraglacial conditions (Figure 6a,e). Under the glacial conditions (Figure 6b,f), on the contrary, the permafrost thickness in MIS-6, which has formed under the glacier, was lower (200 m) than in MIS-2 (360 m). Thus, the glacier impact was expressed in the reduction of the permafrost thickness by 160 m.
During the warm intervals, the permafrost thickness under the extraglacial conditions was lower. As a response to the MIS-5e optimum, it was 210 m. After the warmest interval in the MIS-3 interstadial, because of the sea level rise, it was 270 m at 40 m isobath. At 5 m isobath, warming did not show up at all, since the sea reached only the modern isobath of 40 m (Figure 6e).

Current State of Gas Hydrate Stability Zone
The GHSZ formed together with the permafrost under extraglacial and glacial conditions with heat flux from 50 to 75 mW/m 2 (Figures 9 and 10). The position of the upper boundary of the GHSZ slightly depends on the presence or absence of the glacier in the past. It is located inside permafrost mass a depth of 130-160 m within the extraglacial conditions and at a depth of 140-160 m in the glacial conditions (Figures 9-11).
The lower boundary of the GHSZ within the glacial conditions is located at depths close to 700 m at 50 mW/m 2 , at 450-500 m at 60 mW/m 2 (Figure 12). At 75 mW/m 2 , per our calculations, the hydrate stability zone has not been kept thus far. Under the extraglacial conditions, its depths are 900-1050 and 770-790 m at 50 and 60 mW/m 2 , respectively. At 75 mW/m 2 , the depths are generally close to 400 m. The deep position of the lower boundary and the high GHSZ thickness, especially on 5 m isobath at 50 mW/m 2 , are caused by high permafrost thickness. Accordingly, the lower permafrost thickness in the glacial area determines the shallower position of the lower boundary and the lower thickness of the GHSZ as compared to the extraglacial conditions.

Evolution of Gas Hydrate Stability Zone
The GHSZ evolution is associated with the dynamics of the P-T conditions, i.e., with permafrost dynamics and the glaciation existence. Evolution displays itself, above all, within the periods of GHSZ existence, as well as in the dynamics of its lower and upper boundaries and thickness, caused by climate and sea level fluctuations. Under the extraglacial conditions, the GHSZ lifetime exceeded the estimated period (Figures 9a and 10a). Within the glacial area, the periods of GHSZ changed on a rotating basis with the periods without GHSZ (Figures 9b and 10b). The periods of existence fell on the intervals of 190-125 and from 70-60 ka with a break for the interval of 125-70 ka.
The configuration of the GHSZ lower and upper boundaries and thickness is closely related to the climate, sea level fluctuations and glacier existence. Under the extraglacial conditions, interglacials (MIS-7, 5e, 1) and short-term warm extremums, which were accompanied by the sea level rise (MIS-5c, 5a, 3), on 40 m isobaths are distinguished by depressions of the GHSZ roof up to 100-120 m (Figure 9a). On 5 m isobaths, such depressions are intrinsic only to interglacials (Figure 10a).
The response of the GHSZ lower boundary to the climate fluctuations within the extraglacial conditions repeats the relief of the permafrost lower boundary (Figures 9a and 10a). However, it has a smoother form than the response of the permafrost bottom. The delay is of the same sequence as in the case of the permafrost.
Within the glacial area, two circumstances are noteworthy. The first is the formation of the GHSZ with a roof coinciding with the glacier bed (Figures 9b and 10b). The glacial load causes this. The second is the relief-formed response of the roof, and especially of the GHSZ lower boundary, to the climate and sea level fluctuations on 40 m isobath (Figure 9b). This is not the case on 5 m isobath.
The information on the depth of the GHSZ lower boundary and thickness under a heat flux of 60 mW/m 2 is presented in Table 1. It should be noted that the depth of the GHSZ bottom and its thickness under the extraglacial conditions is significantly more compared to the glacial conditions. Within the pessium MIS-6 they were 800 and 740 m, within the MIS-2-770, and 730 m, respectively. Contrary to the results of the numerical experiment with a heat flux of 60 mW/m2, the calculation based on a flux of 75 mW/m 2 , shows that the gas hydrate stability zone is formed only during the oceanic regressions periods. At the same time, under the glacial conditions, the GHSZ is absent in the recent period. It was destroyed in the Holocene: 4000 years ago in the area of 5 m isobaths, 5000 years ago at 20 m isobaths, and 6000 years ago at 40 m isobaths. Since the conditions favorable for the GHSZ formation are breached, therefore, gas hydrates are possibly subject to degradation methane accumulates in the subpermafrost horizons of bottom sediments. The subsequent destruction of the permafrost layer under the effect of increased heat flux (75 mW/m 2 and more) in the recent period, is likely to cause an increase of bottom deposits permeability and methane flux into the atmosphere from the shallow-water shelf.
The area of high methane emissions northward of the New Siberia Island in the research paper [23] is associated with the existence of the Great Siberian Polynya. The results of our studies indicate a possible connection of this area with permafrost degradation and GHSZ destabilization within the glaciation area of the Middle Neopleistocene with the increased heat flux.

Discussion and Conclusions
In the model, Romanovskii et al. [9] received the values of the permafrost thickness and the depth of the GHSZ lower boundary within the glacial area which were 1.5 times greater than in our model. The differences are explained by the fact that for this research the authors used not conditional, but real geological cross-sections. The latter contains both extended intervals of marine sedimentation, with no rocks freezing, and numerous brown coal layers with low thermal conductivity. Verification of modeling results is an issue of concern in such studies. The concern arises from the unavailability of the data required for comparison. N. Shakhova and D. Nicolsky suggest verification based on the location of the methane seeps [24]. Comparisons are often made with the permafrost thickness in nearshore boreholes [25]. However, the area of the eastern part of the Anjou Islands and the De Long archipelago does not contain such boreholes, except for those used in this study. In addition to testing the scenario and the geological model, the authors used the modeling procedure to make the results more realistic. The calculations for onshore well sites and 5, 20, and 40 m isobaths sites were carried out within a unified system. Additionally, the onshore boreholes sites are considered as benchmarks. This is a geological benchmark: well c-B and the results of generalizing the geological development of the region [55]. In addition, a geothermal benchmark in the form of temperature measurement results for the well c-B. Thus, the modeling results were constantly monitored with the geocryological data in the process of calculations that were considered as verification of the final model.
For the first time, using mathematical modeling, we studied permafrost distribution and thickness within the glaciation area at the end of the Middle Neopleistocene in the northwest of the East Siberian shelf.
It is the first time for the modeling to be based not on conventional schematized sections, but the data of drilling sediments containing, overlying, and underlying the sheet ice of the New Siberia Island. We tested the paleotemperature scenario by the results of the geothermal observations in the boreholes. This ensures the modeling results to be realistic.
The modeling results show that within the glacial area the permafrost thickness at heat flux of 50 and 60 mW/m 2 is 150-200 m less than within extraglacial ones. At the heat flux of 75 mW/m 2 on 20-40 m isobaths, it leads to the formation of discontinuous or sporadic permafrost.
Under the extraglacial conditions, the GHSZ formed at heat flux from 50 to 75 mW/m 2 . Its lower boundary is currently located at a depth of 900-1000 m at 50 mW/m 2 , at a depth of about 400 m at 75 mW/m 2 .
The permafrost and the GHSZ strata continuously exist under the extraglacial conditions, at least during the entire estimated period (200 ka-present). Within the glacial area, during the Middle Pleistocene-Holocene, they occurred and degraded. The recent permafrost strata, as well as GHSZ with the heat flux of 50 and 60 mW/m 2 , exist since MIS-4 (60,000-70,000 years). At the heat flux of 75 mW/m 2 , permafrost is retained only in the coastal area and on 5 m isobaths, and the GHSZ is currently absent.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Physical and thermophysical properties adopted in the numerical model.