Past climate changes and permafrost depth at the Lake El’gygytgyn site: implications from data and thermal modeling

This study focuses on the temperature field observed in boreholes drilled as part of interdisciplinary scientific campaign targeting the El’gygytgyn Crater Lake in NE Russia. Temperature data are available from two sites: the lake borehole 5011-1 located near the center of the lake reaching 400 m depth, and the land borehole 5011-3 at the rim of the lake, with a depth of 140 m. Constraints on permafrost depth and past climate changes are derived from numerical simulation of the thermal regime associated with the lake-related talik structure. The thermal properties of the subsurface needed for these simulations are based on laboratory measurements of representative cores from the quaternary sediments and the underlying impact-affected rock, complemented by further information from geophysical logs and data from published literature. The temperature observations in the lake borehole 5011-1 are dominated by thermal perturbations related to the drilling process, and thus only give reliable values for the lowermost value in the borehole. Undisturbed temperature data recorded over more than two years are available in the 140 m deep land-based borehole 5011-3. The analysis of these observations allows determination of not only the recent mean annual ground surface temperature, but also the ground surface temperature history, though with large uncertainties. Although the depth of this borehole is by far too insufficient for a complete reconstruction of past temperatures back to the Last Glacial Maximum, it still affects the thermal regime, and thus permafrost depth. This effect is constrained by numerical modeling: assuming that the lake borehole observations are hardly influenced by the past changes in surface air temperature, an estimate of steady-state conditions is possible, leading to a meaningful value of 14 ± 5 K for the post-glacial warming. The strong curvature of the temperature data in shallower depths around 60 m can be explained by a comparatively large amplitude of the Little Ice Age (up to 4 K), with low temperatures prevailing far into the 20th century. Other mechanisms, like varying porosity, may also have an influence on the temperature profile, however, our modeling studies imply a major contribution from recent climate changes.


Introduction
The crater Lake El'gygytgyn in NE Russia (67 • 30 N, 172 • 5 E, 492 m above sea level) was formed by an Asteroid impact 3.6 Myr ago (Fig. 1).Since it is believed to have never been covered by ice during the glacial cycles since then, it provides the unique possibility to investigate Arctic climate and environmental changes back to the impact.For this reason, in 2008/2009 an interdisciplinary drilling campaign targeted this area.General information on this International Continental Drilling Program (ICDP) project is given by Melles et al. (2011), while the most important paleoclimatological results were published by Melles et al. (2012).
In the course of the project, boreholes were drilled at two sites, one near the deepest part of the lake (ICDP site 5011-1), and the other on land close to its shoreline (ICDP site 5011-3).From borehole 5011-1, 315 m of sediment cores were retrieved, complemented by 200 m from the underlying volcanic bedrock, while from 5011-3 a 141.5-m-long permafrost core composed of frozen deposits was recovered (Fig. 2).
In this paper, we focus on the characterization of the thermal field beneath and around the lake, studying the influence of variations of thermal properties and past ground surface temperature changes using numerical modeling techniques (e.g.Osterkamp and Gosink, 1991;Galushkin, 1997;Mottaghy and Rath, 2006;Holmén et al., 2011).In the recent years, the impact of climate change on permafrost formation and evolution has become a particular subject of interest.The understanding of the response of permafrost regimes to transient changes in surface temperatures is a key issue regarding the prediction of the influence of global warming on permafrost areas (amongst many others, Saito et al., 2007; Lemke et al., 2007;Serreze et al., 2007;Miller et al., 2010a,b;Serreze and Barry, 2011).
When studying heat transfer in the Earth's upper crust, the upper boundary condition of heat transport is determined by the local climatic conditions.The variations of this boundary condition induces a transient signal which diffuses into the subsurface.Thus, ground temperatures may be seen as an archive of past climate signals.Reconstructing past changes is of major interest, since one of the most important components of climatic change is the variation of temperature at the Earth's surface.However, the diffusive character implies that the older the signal is, the more it is attenuated, with a corresponding larger uncertainty in magnitude and timing.Analyzing the variation of temperature with depth, past fluctuations at the Earth's surface can be reconstructed to a certain extent (see, e.g. the recent review by González-Rouco et al., 2009, and the references therein).This reconstruction can be performed by numerical forward modeling with models of varying complexity, or inverse methods in the narrow sense (amongst others, Shen and Beck, 1991;Beltrami and Mareschal, 1991;Beck et al., 1992;Rath and Mottaghy, 2007).A prerequisite for such a reconstruction is the availability of sufficient data, the most relevant of which is temperature observations.Here, we limit ourselves to forward modeling for determining the present-day temperature profile, because the available temperature data is not sufficient for a reliable reconstruction by inversion of the ground surface temperature history (GSTH) at the study area.Nevertheless, we propose that the combination of forward modeling and sensitivity studies allow some meaningful statements to be made concerning the local GSTH.
As our approach relies on forward modeling, we must assume a certain GSTH entering the simulations.In particular, it has been shown that the amplitude of the Last Glacial Maximum (LGM) and the following warming from Pleistocene to Holocene influence the temperature distribution at all depths up to the surface (Rath et al., 2012).In the arctic region the values of this amplitude are not well known, but the available data suggest a 18 ± 7 K cooler mean temperature than today (Miller et al., 2010a,b).There have been some efforts to determine the spatial distribution of this parameter (Demezhko et al., 2007) whose results also enter our models, but the available data in the area is sparse.In order to do justice to these uncertainties, we perform Monte Carlo (MC) simulations with an ensemble of GSTHs as boundary conditions for the transient forward modeling.At middle-and high-latitude sites like Lake El'gygytgyn, freezing and thawing processes may have a considerable effect on subsurface temperature-depth profiles.Depending on the amount of pores and fractures in rocks, and the fluids filling them, freezing and thawing will play an important role in controlling the subsurface thermal conditions, when cli- matic changes affect the ground surface temperature (Mottaghy and Rath, 2006;Rath and Mottaghy, 2007).This implies that the reconstruction of past temperatures by forward or inverse methods must include these processes.

Available data
The data used for this study originate from the aforementioned boreholes, one in the deepest part of the lake (ICDP site 5011-1, depth 517 m below lake bottom), and one close to its shoreline (ICDP site 5011-3, depth 140 m below surface).Their locations are shown in Fig. 2, which provides a schematic sketch of the basin.Our investigations rely on (a) temperature measurements in both boreholes; and (b) determination of thermal properties from laboratory measurements on cores and γ -ray logging in borehole 5011-1.The temperature data from the lake borehole ICDP site 5011-1 do not represent steady-state conditions.Only the value from the lowermost part is supposed to approximate the undisturbed temperature of the surrounding bedrock (see Sect. 4).Temperature measurements in the permafrost hole 5011-3 started shortly after drilling in November-December 2008.The earliest recordings were expected to reflect conditions undisturbed by the physical effects of the drilling process.However, temperatures around −6 • C were determined at depths that are not affected by daily or seasonal variations.Subsequently, during the years from 2008 to 2010, temperatures were continuously recorded by an automatic data logger in this borehole, in order to monitor the approach of the subsurface temperatures to conditions in equilibrium with the ambient rock.These measurements confirm the first measurements, showing that the value obtained earlier was already rather near to the equilibrium value.In Fig. 3, a selection of the recorded temperatures in the year 2010 is plotted.
It can be seen that at depths of between 15 and 20 m annual variations disappear.Our modeling studies target the impact of long-term climate changes on the temperature field, and consequently annual variations in the uppermost layers are not taken into account.At a depth of 20 m, the average temperature for this year was estimated to be −5.9 ± 0.1 • C.This estimate was in turn used to extrapolate a value for the upper boundary condition of the modeling studies of −6.7 • C at the surface.
Regarding the thermal properties of the prevailing sediments and rocks in frozen condition, we performed thermal conductivity measurements on one core sample retrieved from borehole .
This sample consists of sandy gravel, and is representative for much of the material from the permafrost core (Schwamborn et al., 2012) associated with the model's unit 2 in Fig. 2. We used the full-space and half-space line source measuring device (TK-04, TeKa, Berlin), and the measurement procedures as described by Erbas (2001).Other studies (e.g.Putkonen, 2003) show that the determination of frozen soil thermal properties by a needle probe yields reasonable results.The work flow, elaborated together with the manufacturer, was as follows: first, a part of the frozen core, which consists of loose sand and gravel as well as the ice fraction, was cut off and put in a jar where it was left to thaw.Then, thermal conductivity was determined with the needle probe in the common procedure.To ensure reliable results, several heating cycles were applied.This allowed stable solutions to be found, and also made it possible to determine a mean value along with a standard deviation.
After finishing the measurements on the thawed material, the mixture was frozen with the needle probe remaining in it.Then, thermal conductivity was determined for the frozen mixture.During this procedure, the temperature of the material was ensured to remain well below zero in order to avoid biases due to partial melting around the needle.Finally, the mixture was unfrozen once more, and put in the half-space line source.In contrast to the full-space line source there is less triggering of possible convection since the heating source lies on top of the sample mixture.The results shown in Table 1 confirm this assumed bias: whereas the "unfrozen" values of the full-space line source feature quite some variation, as seen by the standard variation, the values of the frozen mixture are more stable.Furthermore, the thermal conductivity measured by the half-space probe is slightly lower.Porosity was determined by using a defined volume and weighing the saturated and dry sample.This allows determination of the matrix thermal conductivity using the geometric mixing law, which proved to be a good theoretical description of the effective thermal conductivity of rocks with different components (e.g., Hartmann et al., 2008).
Further measurements were performed in order to determine the thermal conductivity of the impact rocks of unit 4 (polymictic and suevitic breccia) and unit 5 (monomictic breccia and fractured bedrock).Two samples were chosen, originating from 334 m (No.C104) and 491 m (No.C170) below lake bottom, respectively.We used the optical scanning method (TCS), which produces a high resolution profile of thermal conductivity along the scanning line of the sample.A detailed description and comparison with other methods can be found in Popov et al. (1999).Thermal conductivity was determined in the original dry state of the sample, after being completely dried, and after being saturated again.The latter two measurement procedures allowed for the estimation of the porosity of the sample.Table 1.Summary of the different measured thermal conductivity values λ (W m −1 K −1 ) of the sample No. 28(2) from the permafrost borehole 5011-3.φ is porosity, which can be estimated from the ratio of volume and mass of the sample (φ a ), or from thermal conductivity (φ b ).Based on these, the corresponding matrix thermal conductivities λ a,b m can be determined.FSLS: full-space line source, HSLS: half-space line source.

Numerical model
The numerical modeling code SHEMAT SUITE (Rath et al., 2006) was used for numerical simulation of heat transport processes at the Lake El'gygytgyn site.As mentioned above, in addition to assuming all thermal properties to be functions of temperatures, latent heat effects due to freezing and thawing of subsurface fluids are accounted for in the code.In particular, the freezing range can be adapted, taking into ac-count that the phase change occurs over a temperature range in the subsurface.Here we assume a functional relationship between unfrozen water content and temperature given by Lunardini (1981) and Lunardini (1988).Water and ice properties in the freezing domain were calculated using the formulations of Fukusako (1990) and Ling and Zhang (2004).
A detailed description of this approach is found in Mottaghy and Rath (2006).3, whereas the magnitude of the Little Ice Age (LIA) and the recent temperature changes are responsible for the curvature in the shallow domain.The mean ground surface temperature before the LIA is taken as a reference.PAT is the preglacial average temperature, TGI is the time of glacial inception, and GAT the average glacial temperature.

Conceptual model
For building the 2-D numerical model, the basin scheme shown in Fig. 2 is used as a basis.This scheme is based on seismic basin profiling (Gebhardt et al., 2006), and coring results (Melles et al., 2011).In Appendix A we give some reasons why the assumption of 2-D geometry is reasonable in the case considered here.
The thermal properties of the different zones were derived from our own measurements (see Sect. 2), as well as literature data for the impact-affected rocks (Popov et al., 2003;Mayr et al., 2008).Core sample 28( 2) is assumed to be representative for the lacustrine sediments of Unit 2, shown in Fig. 2. We deduced some information from the available studies such as Asikainen et al. (2007) on grain size, clay mineralogy, and crystallinity to assess the thermal properties of the lake sediments.Together with a γ -ray log from borehole ICDP site 5011-1 (Gebhardt et al., 2013), and using literature data we estimate thermal conductivity as summarized together with other properties in Table 1.The lacustrine sediments are likely to exhibit lower values of bulk thermal conductivity, which is mainly due to high porosities around 35 % (Melles, 2006).The impact rocks in unit 3 and unit 4 show higher γ -ray values (around 150 gAPI), which implies some quartz content of the suevite bedrock.This is confirmed by the measurements of the two samples from the impact bedrock, and is consistent with results obtained for comparable rocks (Popov et al., 2003;Mayr et al., 2008).Bulk thermal conductivity becomes higher along with lower porosity in the impact rocks.For the transient simulations, additionally volumetric heat capacity must be considered.Here, we use the general mean value of 2.3 MJ m −3 K −1 as given in Beck (1988).With respect to internal heat generation of the rocks, we assume a general value of 0.5 µW m −3 for the lake sediments and 1 µW m −3 for the impact rock, based on the γ -ray data available, as well as literature (e.g.Beardsmore and Cull, 2001).However, this parameter plays only a secondary role since our model only considers the upper few hundred meters of sediments and bedrock, which is not enough for a significant influence of heat generation on the temperature field.Table 3 summarizes some properties of the numerical model.
All thermal properties, in particular the thermal conductivity, represent initial values for the subsequent modeling studies.Only in the deeper part with the impact rocks did Table 4.The five lithological units and their thermal parameters.Thermal conductivity λ is varied within the given ranges, which is ± 20 % of the mean value.φ is porosity, λ bulk denotes the total (bulk) thermal conductivity for the unfrozen material, and ρ c P is volumetric heat capacity.The values in bold correspond to the measured properties of the sample from unit 1, 3, and 4.
the values have to be slightly adapted, using the steady-state simulations presented in Sect. 4 and fitting those results to the temperature measurements in the bottom of borehole 5011-1.

Boundary conditions
The current mean annual air temperature (MAAT) in the area is −10.3 • C (Nolan and Brigham-Grette, 2007;Nolan, 2012), which is a lower bound for the mean annual ground surface temperature (GST).The coupling between air and ground surface temperature is not trivial and has been studied for some decades now (e.g.Smerdon et al., 2004Smerdon et al., , 2006;;Stieglitz and Smerdon, 2007).In particular in areas with varying snow cover, several parameters play a role (Kukkonen, 1987).Generally, the mean annual GST is some degrees above the mean annual air temperature.A recent study (Ge et al., 2011) presents numerical modeling of the different coupling mechanisms between air, soil, and groundwater flow in a permafrost area, the northern central Tibet Plateau.Although PAT is the long term average pre-glacial ground surface temperature, while GAT denotes the minimum temperature of the simplified glacial GSTH.This boxcar-like period of low temperatures begins at TGI.All values here are given as differences, taking the mean ground surface temperature before the LIA as a reference (see Fig. 7).All parameters were assumed to be independently and normally distributed.

PAT GAT TGI T (K)
T (K) kyr this site is not directly comparable to the Lake El'gygytgyn area, the available data confirms a 3-4 K higher mean annual GST compared to the MAAT: as explained in Sect.2, the long-term measurements in borehole 5011-3 allow for the determination of the mean annual GST at the study site, which is about −6 • C at 20 m depth.
For the transient simulations we use a GSTH as discussed in Sect. 5.At the lake floor, a mean temperature of 5.5 • C is set, which is supported by some temperature data (see Sect. 3).Whereas the bottom lake water approaches 4 • C, the value where water density is highest, the higher temperature value in the uppermost lake sediment layers may point to a seasonal signal that has been preserved from previously warmed summer water masses.It is known from subglacial lakes in Antarctica without seasonal cycling that water temperatures at the bottom of this kind of lake remain at 4 • C at maximum, with only slight variation (Vincent and Laybourn-Parry, 2008).The water levels of Lake El'gygytgyn have changed on long-term scales.For instance, it is well documented that the level was some 10 m lower at the LGM (Juschus et al., 2011).However, this change in depth is small compared to the total water depth of up to 175 m in Lake El'gygytgyn, and probably did not influence the water mass circulation during the last 125 kyr.Therefore, we accept the measured value of 5.5 • C as a reasonable value for the upper lake sediment environment (see Sect. 4).
The lower boundary condition is an estimated basal heat flow, as there is no data in the vicinity of the study area.The compilation of global heat flow data of Pollack et al. (1993) lists some sparse values originating from boreholes some 500 km away.Yershov (1998), referring to Russian sources, presents a heat flow map with values between 40 mW m −2 and 80 mW m −2 for this area.Recently Goutorbe et al. (2011) published a new compilation of global heat flow, combining the data of Pollack et al. (1993) with geological proxies, applying different interpolation techniques.From their findings, the surface heat flow in the area is 66 ± 8 mW m −2 .Therefore, we initially set the lower boundary condition to this value and compare the results to the available temperature data (see Sect. 4).

Steady-state simulations
In order to determine the temperature field near and beneath the lake, we performed steady-state simulations as a first step.All the parameters described in Sect.3.2 enter the numerical modeling.In particular, both the temperature profile in the permafrost borehole, as well as the lake borehole, give valuable constraints to the upper and lower boundary conditions.However, the data from the lake do not represent undisturbed conditions.Thus, the data is restricted to some general temperature gradients, because the temperature measured in the lowermost part (400 m below lake surface) should be most reliable (Kueck, personal communication, 2010).The different logs and the assumed gradient are shown in Fig. 4. Therefore, we used this value (about 25 • C) to test the initially assumed basal heat flow (see Sect. 3.2).However, it only slightly varied, changing from an initial value of 66 mW m −2 to 70 mW m −2 , which is still within the uncertainty range given in Goutorbe et al. (2011).Therefore we can conclude that the input parameters of our model represent a reasonable basis for the subsequent modeling studies.
In contrast to the lake hole temperature data, the long-term record in borehole 5011-3 reflects undisturbed conditions, but the shallow depth (140 m) only allows a limited calibration with respect to the thermal properties of the layers beneath.Furthermore, transient effects such as paleoclimatic ground surface temperature changes are likely to have a significant influence on the temperature field, as described in Sect. 5. Thus, the measurements of thermal conductivity give valuable constraints for reliable temperature modeling.
Figure 5 shows the modeled temperature profiles in comparison with the data plotted in Fig. 4 at borehole 5011-1 in the lake, as well as at the land location, borehole 5011-3.The dashed and dotted lines show the strong sensitivity to thermal conductivity.The values are varied within the ranges given in Table 4.The lower and higher thermal conductivity values yield a higher and lower temperature gradient, respectively.This variation in thermal conductivity is rather large: due to  10.The mean ground surface temperature before the LIA is taken as a reference (see Fig. 7).6. Top: Line plots of the calculated subsurface temperatures, with the median highlighted.Bottom: corresponding density plot and average temperature profile for the same data set.
our measurements, the actual values are likely to be close the mean values, represented by the solid lines in Fig. 5.
It is important to note that steady-state simulations can explain the temperature distribution beneath the lake, applying thermal conductivity values determined as described in Sect.3.1.As explained in Sect.3.2, the lowermost temperature value is considered to be representative and can be well reproduced by the model.However, the steady-state simulations are not able to explain the temperature profile of the permafrost borehole 5011-3 (Fig. 5).Even a variation within a rather large uncertainty range of ±20 % is insufficient to match the measured data.Therefore, transient surface temperature changes must be taken into account.
In  10.This depth is defined at a temperature of 0 • C, which is a lower limit for this parameter.
temperature field to this parameter.The lower boundary of permafrost (0 • C isoline) is highlighted in order to visualize the sensitivity to thermal conductivity.

Transient simulations
In order to study the influence of past climate variations, we performed transient simulations, applying a simple boxcar model representing climate changes over the past 125 kyr (Fig. 7).For this, we used the available information about the past climate as described in Sect. 1. From pollen records, only vegetation and thus information about summer climate conditions can be reconstructed, which is obviously only a partial contribution to the mean annual air temperature.Furthermore, the duration, thickness and kind of snow cover play a role in the coupling between air and ground surface temperatures (see Sect. 3.2).As the study site seems to have remained without glacier cover, very low temperatures during the LGM are likely to have prevailed (e.g., Miller et al., 2010a).Therefore, a rather large amplitude of the Pleistocene-Holocene temperature step is likely.
The results of our transient simulations confirm this hypothesis: in order to explain the temperature variation with depth in the lower part of the measured profile, the amplitude of the LGM (denoted as PGW in Fig. 7) is estimated to be around 14 K, which allows the model to best reproduce the temperature gradient at the base of the observational profile in 5011-3.This agrees well with results from earlier studies (Miller et al., 2010b,a).The values of LGM temperatures derived from the study of Melles et al. (2012) are not directly comparable, as they only refer to maximum summer temperatures.However, the value of 14 K determined from our modeling is considerably higher than those recently determined by other authors (Schneider von Deimling et al., 2006;Schmittner et al., 2011;Annan and Hargreaves, 2012).Yet, some uncertainty remains in our results because of the insufficient depth of the permafrost borehole: the available short temperature profile only allows determination of the temperature gradient, which is influenced by the basal heat flow, thermal conductivity, and past climate changes.Nevertheless, we have some control on basal heat flow and thermal conductivity based on the temperature measurements in the deepest part of borehole 5011-1, as well as by the thermal conductivity measurements.The temperature observations in 5011-1, located far away from the limits of the lake, should not to be influenced by SAT changes due to the overlying water body.Thus it can be used for an independent determination of steady-state heat flow.From this reasoning, we consider the relatively large amplitude of post-LGM warming as probable.In order to account for the ambiguity of heat flow density and past climate changes, we also performed several sensitivity studies, which showed that that our results have an uncertainty of about ±5 K in the amplitude of the LGM.
The results of the transient simulations are shown in Fig. 8.This figure also shows the situation at the end of the LGM when permafrost depth was obviously much larger.Today, our modeling results suggest a permafrost base of approximately 340 m depth at the borehole 5011-3.Compared to the results from the steady-state simulations, it becomes evident that the past changes in GST have a significant influence on the temperature distribution today.Table 5 summarizes the different results for this value at the permafrost borehole 5011-3.
The upper part of the measured profile in borehole 5011-3 shows some strong curvature (see blue lines in Fig. 8), which was thoroughly investigated by sensitivity studies in terms of more recent surface temperature changes.Other possible influences are unknown heterogeneities in thermal parameters, porosity changes, or fluid flow.However, these can largely be excluded due to several reasons: on the one hand thermal properties have been determined by different approaches (indirect as well as measurements), on the other hand permafrost has been present for at least some ten thousand years.The warmer period during the Allerød (13.3-13.1 kyr BP), when the lake level was some 5 m higher than modern levels and flooded the coring site at that time (Schwamborn et al., 2008(Schwamborn et al., , 2012)), was studied by different model runs and showed that there is no influence on the current temperature field anymore.
Therefore we conclude that past surface temperature changes are most likely to be responsible for the temperature variation with depth.As a result we propose a possible large amplitude of the LIA of 3.5 K lower temperature during the 1500-1900 period GST relative to today's mean GST.A subsequent still cool period until 1950 (of 2.5 K lower than today) is able to approximately explain the observed temperature profile.Though data are sparse in this area of Siberia, there are some other indications for low LIA temperatures.Mann et al. (2009) include some tree ring density reconstructions from northern Siberia, which indicate temperature differences of more than −2.5 K with respect to modern temperatures.In addition, the 10-yr averaged land surface temperature anomalies from the recently published CRUTEM4 database (Jones et al., 2012), which include many new data from Siberia, show a strong warming of more than 1 K since 1970 in this area.
In order to obtain an overview of the uncertainties of our estimate, we performed an MC experiment.The setup was as following: as in the simulation presented above, the last glacial cycle was parametrized as simply as possible, in this case a boxcar function.As the time of the postglacial warming is rather well constrained, we choose to vary only three parameters, the preglacial average temperature (PAT), the average glacial temperature (GAT) and the timing of the glacial initiation (TGI).Concerning the PAT, it is well known that due to the strong attenuation of temperature perturbations with time, only the temperature changes of the last glacial cycle play a significant role.Older variations may safely be approximated by an appropriate average temperature (see, e.g.Majorowicz et al., 2008).For the purposes followed here, even the variations during the glacial cycle are accommodated easily this way.Though nearer to the present, even the steepness of the temperature rise during the glacial-interglacial transition plays a detectable but minor role.In particular, Rath et al. (2012) have shown that for very shallow boreholes as the one considered here, the joint effect of the glacial-interglacial complex is indistinguishable from an additional linear trend, or, more physically, to an additional heat flow density (see above).Thus, we think this simple parametrization is reasonable.
We assumed independent normal distributions for all three parameters, implying that they are characterized by their means μi and standard deviations σi .The values assumed are given in Table 6, while the generated distributions can be found in Fig. 9 results we conclude that the uncertainties assumed for the more ancient GSTH produce a considerable spread in temperatures already in the first few hundreds of meters.This spread carries over to the depths of the lower limit of permafrost, which can be found in Fig. 11.Unfortunately the existing geophysical surveys do not extend beyond the lake area and its immediate surroundings, so that recent undisturbed permafrost extent could not be observed.However, the values in Fig. 11 largely agree with the values found in the literature for this area (Yershov, 1998).

Conclusions and outlook
Our results of 2-D numerical simulations of the temperature field in El'gygytgyn Crater indicate the permafrost depth and the talik dimensions around and below the lake basin, as well as the influence of past temperature changes.The temperature model suggests a permafrost depth surrounding the lake of 330-360 m.This result depends on the assumption of an ancient transient signal present in the temperature distribution that evolved from low temperatures during the Last Glacial Maximum.Sensitivity studies and MC simulations permitted a crude characterization of the uncertainties in permafrost depth.More recent signals in the temperature distribution presumably originate from the LIA and the recent global warming.In particular regarding the LIA, our comparison between data and modeling suggests a significant and persistent cooling during that time period.A possible bias by fluid flow can probably be excluded, because the regime below the active zone remained frozen for at least the last 10 kyr.The talik below the lake basin is pervasive and follows a normal tem-perature gradient around 3 K per 100 m, starting at 5.5 • C in the uppermost lake sediment layers.
Future studies are necessary to improve the interpretation of the very shallow temperature field.In particular, these will include additional forward modeling in order to improve our conceptual model, and assess the sensitivities of our results with respect to subsurface heterogeneity or possibly subsurface flow.Using more sophisticated inverse techniques will not only enable us to obtain a better characterization of the past surface temperature changes, but also to characterize the uncertainties involved.
The annual variation visible in the high quality temperature monitoring data of the top 20 m was not investigated in this study.However, these observations could be used to retrieve valuable constraints for the very shallow subsurface properties.

Appendix A Is 2-D modeling adequate?
The assumption of a 2-D geometry for Lake El'gygytgyn and its surroundings seems to put severe constraints on the validity of the results presented here.However, due to the particular geometry found (large lake diameter of ≈ 12 km, small lake depth of ≈ 170 m) this simplification turns out to be justified.In order to assess possible biases introduced by the 2-D approach, we set up a simplified 3D model with a rather realistic shape of the lake, but with a constant average depth of 170 m.Furthermore, it only consists of two geological layers, one with a higher porosity representing the lake sediments and one with a lower porosity representing the impact-affected bedrock.For this simplified setup, we followed the same simulation procedure as described in Sects.4 and 5.In order to be able to compare possible differences between a 2-D and a 3-D approach directly, we additionally simulated a 2-D model crossing the location of the permafrost borehole (see Fig. 12a).The results are shown in Fig. 12b and c.Differences in temperatures between the 2-D slice extracted from the 3-D model and the 2-D model itself are small, with a maximum value of about ±1 K.These extreme values appear beneath the lake boundary, which is to be expected since the largest temperature gradients occur here.We observe only a slight difference in temperature at depths comparable with the borehole, shown in Fig. 12.The small deviation between 2-D and 3-D model results yields an uncertainty of about 10 m in the lower bound of permafrost.Hence, we conclude that the restriction to a 2-D model is justifiable for this investigation.
Fig. 1.(a) Geographic position of the El'gygytgyn Crater in NE Russia (red box).(b) A Landsat image showing the crater (dotted line) and the drill sites 5011-1 and 5011-3.

Fig. 2 .
Fig. 2. Basin scheme of the Lake El'gygytgyn site with the digital elevation model in the background.The numerical model is based on this conceptual model by Melles et al. (2011).

Fig. 3 .
Fig. 3. Selection of the temperatures logs in the permafrost borehole 5011-3 during the year 2010.Note the different scales in the two panels.

Fig. 4 .
Fig. 4. Temperature logs in borehole ICDP site 5011-1.T 5 is the most reliable, and thus is used for estimating the temperature gradient (data by ICDP-OSG, 2009).

Fig. 5 .
Fig. 5. Modelled temperatures at borehole locations 5011-1 and 5011-3.The dashed lines show the sensitivity to thermal conductivity variations.The colours indicate the units as in Fig. 2.

Fig. 6 .
Fig. 6.(a) The model and its units.The temperature field based on the steady-state simulations using (b) the mean value, (c) lower thermal conductivity, and (d) higher thermal conductivity.The vertical scale is greatly exaggerated.

Fig. 7 .
Fig. 7.The simplified boxcar model used in the transient simulations.The magnitude of the post-glacial warming (PGW) is influencing the gradient of the deeper part of the temperature profile in Fig.3, whereas the magnitude of the Little Ice Age (LIA) and the recent temperature changes are responsible for the curvature in the shallow domain.The mean ground surface temperature before the LIA is taken as a reference.PAT is the preglacial average temperature, TGI is the time of glacial inception, and GAT the average glacial temperature.

Fig. 8 .
Fig. 8. Results from the transient simulations.Temperature profiles at borehole 5011-3 showing (a) the whole model depth, and (b) only the upper part.The temperature field of the model (c) today, and 20 kyr BP (d).The vertical scale is greatly exaggerated.

Fig. 9 .
Fig.9.Histograms of the parameter values for the Monte Carlo simulations in Figure10.The mean ground surface temperature before the LIA is taken as a reference (see Fig.7).

Fig. 10 .
Fig. 10.Results from the MC simulations for the 5011-3 borehole.The parameter values assumed are given in Table6.Top: Line plots of the calculated subsurface temperatures, with the median highlighted.Bottom: corresponding density plot and average temperature profile for the same data set.

Fig. 11 .
Fig.11.Histograms of the depth to the lower boundary of permafrost for the Monte Carlo simulations in Fig.10.This depth is defined at a temperature of 0 • C, which is a lower limit for this parameter.

Fig. 12 .
Fig. 12.Comparison of temperatures derived from simplified 3-D and 2-D models for the lake area.(a) Simplified 3D model for Lake El'gygytgyn and its surroundings, based on a quasi-realistic representation of the lake in its horizontal extension, but with a constant depth of 170 m.Also shown is the position of the vertical slice chosen for the 2-D model, and the 5011-3 borehole.(b) Temperature differences (K) between a vertical slice through the full 3-D model, and the corresponding 2-D model.(c) Modeled borehole temperature profiles at the 5011-3 borehole.

Fig. 13 .
Fig. 13.Results from the model runs with varying porosity.The largest influence occurs in the temperature range where latent heat it is released or absorbed.
Table 2 lists the results.

Table 2 .
Summary of the different measured thermal conductivity values λ (W m −1 K −1 ) for the two samples from unit 3 (breccia) and unit 4 (fractured bedrock), respectively, minimum (min) and maximum (max) values along the scanning line, as well as mean values.Porosity φ and matrix thermal conductivity λ matrix are determined from saturated λ sat and dry thermal conductivity λ dry measurements.λdry λ min-max,dry λ sat λ min-max,sat φ from λ (-) λ matrix

Table 3 .
Properties of the 2-D model.In case of the transient modeling, the temperature at the upper boundary is time dependent as shown in Figure7.

Table 6 .
Parameters used for the Monte Carlo simulations in this study.