Thermal conductivity of firn at Lomonosovfonna, Svalbard, derived from subsurface temperature measurements
Accurate description of snow and firn processes is necessary for estimating the fraction of glacier surface melt that contributes to runoff. Most processes in snow and firn are to a great extent controlled by the temperature therein and in the absence of liquid water, the temperature evolution is dominated by the conductive heat exchange. The latter is controlled by the effective thermal conductivity k. Here we reconstruct the effective thermal conductivity of firn at Lomonosovfonna, Svalbard, using an optimization routine minimizing the misfit between simulated and measured subsurface temperatures and densities. The optimized k* values in the range from 0.2 to 1.6 W (m K)−1 increase downwards and over time. The results are supported by uncertainty quantification experiments, according to which k* is most sensitive to systematic errors in empirical temperature values and their estimated depths, particularly in the lower part of the vertical profile. Compared to commonly used density-based parameterizations, our k values are consistently larger, suggesting a faster conductive heat exchange in firn.
Glaciers and ice sheets are important indicators of past and ongoing climate changes. Under the influence of temperature fluctuations at the surface the subsurface glacier temperature also changes. As a basic physical property of a medium, temperature of snow, firn and ice controls multiple processes occurring therein and at the glacier surface. Climate-induced glacier mass change is strongly affected by the state of snow and firn, where liquid water generated at the surface of glaciers during the ablation period can be refrozen, thus reducing runoff. The magnitude of liquid water refreezing is largely dependent on the subsurface temperature (e.g. Trabant and Mayo, 1985). A snowpack or firn pack reaching lower temperatures during the winter season is able to refreeze a larger amount of water during and after the ablation season. Warmer snow and firn experience faster metamorphism (e.g. Jordan et al., 2008) and gravitational densification (Ligtenberg et al., 2011). The distribution of temperature under the glacier surface also defines the ground heat flux, which contributes to the energy balance at the surface, and is thus important for simulation of the surface energy fluxes and melt rates. In addition, cold ice is more viscous and less prone to deformation (e.g. Weertman, 1983). Therefore, a colder glacier can be expected to exhibit lower flow velocities, provided that other environmental parameters are equal. Thus the processes of mass and energy exchange occurring at glaciers are in tight interaction with the subsurface temperature, which needs to be either measured or simulated, if measurements are not possible.
In the absence of liquid water the changes in subsurface temperature are defined by the process of thermal conduction described by Fourier's law (Cuffey and Paterson, 2010, p. 403), according to which the heat flux (Q) is proportional to its thermal conductivity (k) and to the spatial temperature gradient (∇T): . Since most temperature fluctuations occur at the surface, the dominant direction of heat flux in near-surface glacier layers is vertical: , where z is the vertical coordinate. Sturm et al. (1997) indicated that several processes contribute to the temperature changes occurring in a subfreezing porous snow and firn, namely conduction through the rigid ice matrix, conduction through the air in pores, and latent heat transport through the pores due to sublimation and condensation of water vapour. To underline that fact and because all three processes are essentially driven by the temperature gradient (e.g. Bartelt and Lehning, 2002), here we use the term effective thermal conductivity (k) following Sturm et al. (1997) to describe the ability of snow, firn and ice to transport thermal energy. Along with density (ρ) and specific heat capacity (C) k is used to calculate thermal diffusivity (κ) as .
Due to the direct connection between k of a medium and temperature changes therein, most empirical estimates of k are based on temperature measurements (Sturm et al., 1997). Continuous measurements of natural temperature fluctuations in snow and firn allow us to derive κ and k values using either Fourier-type analysis or optimization techniques. Diurnal temperature fluctuations penetrate down to ca. 1 m and the associated phase lag and/or amplitude dampening occurring with depth can be used to reconstruct the effective thermal conductivity of seasonal snow (e.g. Sergienko et al., 2008; Osokin and Sosnovsky, 2014). Annual fluctuations penetrate down to ca. 10 m and thus k estimates for thick firn packs using the method require a long data series undisturbed by the influence of liquid water (Dalrymple et al., 1966; Weller and Schwerdtfeger, 1971; Giese and Hawley, 2015). Alternatively, κ and k values in the heat equation can be determined by minimizing the misfit between simulated and measured natural evolution of temperature in snow and firn (e.g. Tervola, 1989; Zhang and Osterkamp, 1995; Yang, 1998). Brandt and Warren (1997) applied this method for near-surface snow at the South Pole, Sergienko et al. (2008) for the snowpack on a drifting iceberg and Oldroyd et al. (2013) for the seasonal snowpack on an Alpine glacier.
Another option is to induce a heat flux in the snowpack using a heat source and register the temperature response in either an object with well-known properties that is in contact with both the heat source and the snow or snowpack itself. The former method is known as the needle probe technique and is widely used to measure in situ effective thermal conductivity of porous materials including snow (e.g. Lange, 1985; Singh, 1999; Morin et al., 2010). In the latter case the snow sample is placed on a heated plate controlling the vertical heat flux and the respective temperature gradient is measured, and the relation between the two values yields k (e.g. Calonne et al., 2011; Riche and Schneebeli, 2013). Sturm et al. (1997) provided an extensive overview of the above-named methods and associated uncertainties, which was followed up by further insights into possible biases of needle probe measurements (Riche and Schneebeli, 2010).
A novel technique for estimation of k of snow was suggested by Kaempfer et al. (2005) and further developed by Calonne et al. (2011) and Riche and Schneebeli (2013). The method relies on numerical simulation of processes contributing to the heat flux driven by the temperature gradient based on detailed 3-D X-ray micro-tomography images of the snow matrix. It allows us to obtain the k tensor for relatively small snow samples.
Published k values for snow and firn vary from 0.1 to 2.5 W (m K)−1 with a strong correlation with density, which justifies the use of density as a proxy to calculate the effective thermal conductivity for modelling purposes (e.g. Sturm et al., 1997; Riche and Schneebeli, 2013). The considerable spread in values suggested by different parameterizations is explained by the inconsistency in applied measurement techniques and associated uncertainties and also by the influence of snow and firn parameters other than density. That can be temperature, grain size and contact area, pore diameter, and inter-connectivity and anisotropy of the k property.
The purpose of the present study is to reconstruct the values of effective thermal conductivity of a thick snowpack and firn pack at Lomonosovfonna, Svalbard, based on evolution of subsurface temperature and firn density measured in 2012–2015. Effective thermal conductivity of firn is derived for five distinct periods by minimizing the misfit between the measured and simulated subsurface temperature evolution. The method is promising, particularly for thick firn packs. To our knowledge it has not been applied for this purpose so far, although Sergienko et al. (2008) employed similar routines for a seasonal snowpack and Nicolsky et al. (2007) for permafrost. Fourier analysis applied earlier for a thick firn pack at the Summit Station of the Greenland Ice Sheet (Giese and Hawley, 2015) cannot be used here due to the influence of meltwater. Retrieval of snow and deep firn samples for direct measurements using heated plate or needle probe methods is logistically challenging. Estimates of the firn k values are complemented by uncertainty quantification experiments exploring propagation of possible biases in empirical data through the applied models.
We use data on subsurface density and temperature evolution collected at Lomonosovfonna, Svalbard, a flat ice field nourishing several outlet glaciers. The field site is at 78.824∘ N, 17.432∘ E, 1200 m a.s.l., which is well above the equilibrium line, estimated to be at ca. 720 m a.s.l. (van Pelt et al., 2012). The local glacier thickness is 192±5.1 m (van Pelt et al., 2013) of which the firn layer constitutes ca. 20 m (Wendl, 2014). The accumulation rates estimated from repeated radar surveys (Pälli et al., 2002; van Pelt et al., 2014) are 0.58–0.75 m w.e. yr−1. The melt rate simulated by a model describing surface fluxes of energy and mass (van Pelt et al., 2012) is 0.34 m w.e. yr−1. The net annual accumulation rates resulting in the relative vertical shifts of the glacier surface are estimated at 1.12 (April 2012–2013), 1.32 (April 2013–2014) and 0.9 (April 2014–2015) m. The first two values are based on readings at mass balance stake S11 (see Fig. 1a in Marchenko et al., 2017a) and the last one is derived by minimizing the misfit between temperature profiles measured in April 2015 by the thermistor strings installed in April 2014 and in April 2015. The firn pack at Lomonosovfonna is heavily influenced by the percolation and refreezing of meltwater (Marchenko et al., 2017b), which results in prominent variability of subsurface stratigraphy at the scale of 10 m (Marchenko et al., 2017a).
2.1 Subsurface density
Four shallow firn cores (9.5 cm diameter) were drilled at Lomonosovfonna in April 2012–2015 using a Kovacs core drill. The density and stratigraphy profiles (Fig. 1) for the first three cores were presented in detail by Marchenko et al. (2017a) along with the details on field and laboratory procedures applied. Similar routines were applied in 2015.
2.2 Subsurface temperature
Subsurface temperature was measured using multiple thermistors grouped in several strings. They were placed in boreholes with a 5.5 cm diameter drilled by a Kovacs auger. The boreholes were then backfilled with drill chips and loose snow to minimize the perturbation of the snow and firn media. In 2012–2014 nine thermistor strings were placed in a rectangular 3×3 grid with a separation of 3 m between neighbouring boreholes (Marchenko et al., 2017a, b). In April 2015 only one thermistor string was installed.
All thermistor strings were custom manufactured at Uppsala University using a multi-leaded cable with a PVC jacket and uni-curve NTC thermistors. In 2015 the cable was placed in multiple 2 m long rigid plastic tubes to ensure a precise and constant separation between neighbouring sensors. The sensors were fed through holes in the tubes and fixed at their outer surface for better contact with the sounded environment.
The thermistor strings were scanned using a Campbell Scientific CR10X data logger and several relay multiplexors. Each thermistor was connected in series with a reference resistor and precisely measured excitation voltage was applied to the circuit. Voltage drop over the reference resistor was measured and then converted first to corresponding resistances of the sensors and then to temperature values using Ohm's law and recommendations provided by the thermistor manufacturer.
The evolution of subsurface temperature was recorded during four time periods: 21 April–19 October 2012, 22 April–12 July 2013, 17 April 2014–11 April 2015 and 15 April–9 July 2015. The frequency of measurements is shown by the colour bar at the bottom of Fig. 2. During the first two periods it was once every 3 h; during the fourth period it was once every 1 h. During the third period the frequency varied and was once every 1 h during 17 April–31 July 2014 and 15 April–9 July 2015, once every 3 h during 1 August–31 October 2014 and once every 12 h from 1 November 2014 to 14 April 2015. The strings installed in 2012, 2013 and 2014 contained up to 128 thermistors covering a depth from 0.5 to 12 m with a vertical separation varying from 0.25 to 2 m (Marchenko et al., 2017b; see Fig. 2 therein). In 2015 the string contained 31 sensors separated by 0.25–1 m covering a similar depth interval.
2.2.2 Data post-processing
Several post-processing routines were applied to the measurements of the subsurface temperature evolution. Firstly, the individual thermistors were calibrated against 0 ∘C by applying an offset defined as the mode (most frequent value) of the values measured during 1 July–1 September. During this period subsurface temperature is most likely to reach 0 ∘C. The value serves as a natural upper limit for the snow and firn temperature and is well interpreted from plots of measured values against time.
Secondly, spline interpolation was applied to interpolate the data from each thermistor string to a common vertical grid with 0.1 m spacing between neighbouring nodes. During April 2012–2015 the subsurface temperature was simultaneously measured by several thermistor strings and thus the interpolated dataset was subsequently averaged horizontally to produce a single series τ describing subsurface temperature changes in time (superscript n) and depth (subscript i):
where is the temperature measured by the mth thermistor string and q is the total number of strings. In 2015 only one thermistor string was installed and thus no averaging was applied. To characterize the spread in temperature values measured at the same depth and time but by different instruments the corresponding standard deviations were calculated as
The evolution of subsurface temperature in the upper 10–15 m of a glacier is mostly controlled by two processes: conductive heat flow and refreezing of liquid water accompanied by release of the latent heat (Cuffey and Paterson, 2010, p. 403). Configuring the computational domains to minimize the influence of non-conductive heat fluxes, we compute the effective thermal conductivities for the firn profile at Lomonosovfonna and assess their sensitivity to errors in empirical data used in the simulations.
The computational procedure is as follows. Within the forward model (see Sect. 3.1) the heat equation is approximated numerically and then solved for the temperature with a given conductivity k and density ρ. The two parameters are then iteratively adjusted to derive conductivity k* and density ρ* minimizing the difference between the simulated and measured subsurface temperature and density (see Sect. 3.2). To quantify the uncertainties associated with k* we first define the feedback of simulated temperature to change in individual conductivity values. These results are then used to compute the sensitivity of conductivity to errors in the empirically derived
temperature values (see Sect. 3.3.1),
depths of temperature values (see Sect. 3.3.2) and
density values (see Sect. 3.3.3).
These routines were coded using MATLAB R-2018a software. All computations on a laptop with four cores Intel(R) i7-4600 CPU 2.1 GHz are obtained in less than 10 s.
The time periods are chosen to minimize the influence of water refreezing on the evolution of subsurface temperature. They are referred to as “spring 2012”, “spring 2013”, “spring 2014”, “autumn 2014”, and “spring 2015” and correspondingly cover 21 April–19 June 2012, 22 April–1 June 2013, 18 April–4 July 2014, 25 September 2014–11 April 2015 and 15 April–9 July 2015. Furthermore, temperatures values above −2 ∘C are excluded from the analysis to avoid the influence of latent heat fluxes from firn volumes with increased water content as the freezing front propagates through them. To minimize the potential influence of near-surface processes (radiation penetration, wind effects) on the results, we disregard the temperature values measured above the depth of 1 m referenced to the glacier surface at the moment of instrument installation.
3.1 Conduction model setup
The model is based on Fourier's law of heat conduction. The temperature of the firn T (∘C) at depth z (m) and time t (s) is governed by the one-dimensional equation
where ρ is the subsurface density (kg m−3), C=2027 J (kg K)−1 is the specific heat capacity calculated using the temperature-dependent function from Cuffey and Paterson (2010, p. 400) for the temperature of −10 ∘C and k is the effective thermal conductivity (W (m K)−1). Given k(z), ρ(z), and the initial and boundary conditions, Eq. (3) is solved forward in time for T(z, t).
The numerical solution of Eq. (3) is based on a discretization with the time step Δt and space step Δz. The spatial and temporal derivatives are approximated using the Crank–Nicolson method (Dahlquist and Björck, 2003) with central finite differences in space and trapezoidal rule in time. Let be the temperature T(zi, tn) at tn=nΔt, and zi=iΔz, . The number of time steps is N and the number of nodes in space Mn varies in time. The solution at zi is advanced in time from tn to tn+1 as
where is the spatial temperature derivative at zi approximated by
Collecting terms in Eq. (4) with Tn+1 on the left-hand side and terms with Tn on the right-hand side, a tridiagonal system of linear equations can de derived and solved for by a direct method. The time integration is unconditionally stable and second-order accurate in space and time (Dahlquist and Björck, 2003).
At each time step the computations are performed for Mn nodes from the top at z1=1 m to the bottom at where the temperature is just below −2 ∘C, which minimizes the influence of the latent heat release at the freezing front. Thus the depth of the lower boundaries of the simulation domains increases over time following the downward propagation of the −2 ∘C isotherm as shown by the white curves in Fig. 2. The vertical step (Δz) follows the depth interval in the empirical dataset and is 0.1 m. The chosen time step is 1 h and in case the measurement period is larger than that, linear interpolation is used to derive the values missing in the upper and lower boundary conditions. The model is initialized in each of the five time periods using measured temperature values . The upper and lower boundary conditions are of Dirichlet type and are determined by the temperatures measured at z1 and : and .
The density ρ(z) and effective thermal conductivity k(z) at depths zi are constrained using piecewise linear interpolation based on J nodes vertically spaced by 1 m. Since the forward model is used within an inversion routine (see Sect. 3.2) optimizing the J 1 m spaced ρ and k values, the latter are given by an arbitrary initial guess at the first inversion iteration and later by the results of the previous inversion iteration. For the domains covering spring seasons J equals 8 and is 6 in autumn 2014 with the uppermost value corresponding to the upper boundary of the computational domain. The choice of J value is a compromise between coarse vertical resolution of the optimized parameters (low J) and insufficient number of data to constrain a very detailed k* profile (high J). Too large J will result in an oscillatory optimal solution for k.
3.2 Inverse routine
The effective thermal conductivity k in Eq. (5) is unknown. Therefore, the above-described forward model is used in an optimization routine to iteratively derive the values of the effective thermal conductivity and density that minimize objective function . Following Smith (2013), the latter is defined by the sums of squared differences between the simulated and measured temperature and density:
and the optimization routine attempts to minimize by adjusting the conductivity k and density ρ values used in simulations. In Eq. (6) and are the simulated and measured temperature values at time tn and depth zi, and ϱi is the density measured at depth zi. For the spring domains ϱ profiles are taken from the cores drilled within a couple of days from the start of simulation and for the autumn 2014 domain we reuse the density data from the core drilled in April 2014. The deviations of the simulated temperature values from empirical data are weighted by the variances in temperature values from different thermistor strings but at the same time and depth (see Eq. 2). This results in lower significance of simulation errors when the measurements are less certain. It is thus assumed that measurement errors are independent and normally distributed with zero mean and variance . Since only one shallow core was drilled every year, empirical data are not available to quantify the errors in density measurements, and the weighting term was set to 1 kg m−3. The value γ=10 is chosen to keep the balance between the temperature and density terms in F such that the optimal solutions are smooth and ρ is close to the measurements in l2-norm. The choice of γ depends on the size of the data N, Mn, and the magnitudes of τ and ϱ.
Choice of the cost function in Eq. (6) does not assume any correlation between k and ρ. This relation is derived later in Sect. 4.6 based on the optimized k and ρ values. The primary aim here is to derive the optimal k values at J nodes. Optimization of the densities used by the forward model is included in the optimization to allow for flexibility in the parameter ρ since the measurements in single firn cores are uncertain and may be not representative at the scale of ca. 10 m covered by thermistor measurements. The second term on the right-hand side of Eq. (6) can be interpreted as a Bayesian prior guess (Smith, 2013; Calvetti and Somersalo, 2007) adding extra information about ρ or a regularization of the density according to Tikhonov (Calvetti and Somersalo, 2007).
We introduce the vectors:
where . Let the diagonal matrices W have , , , , and Wjj=0, , , , on the diagonal. Then F in Eq. (6) can be written as
The diagonal elements of W vanish when i>Mn since the sums in Eq. (6) are restricted to i≤Mn.
It is solved by the lsqnonlin function in MATLAB.
3.3 Uncertainty quantification
Results of the optimization routine largely rely on the empirical data used to guide the routine. The sensitivity of optimized k* values is explored separately for the errors in measured temperature, depths of empirically derived temperature values (both affecting the vertical temperature gradient) and measured density values by applying uncertainty quantification techniques described by Smith (2013). The results of the sensitivity estimates are only valid for relatively small errors on the order of ca. 5–10 % of the parameter value in question.
Errors δτ in the measured temperature values τ in Eq. (6) propagate to the effective thermal conductivities k* derived from the optimization problem in Eq. (8). Here we first inverse the logic and calculate the feedback of simulated temperature values T to relatively small perturbations in individual k* values. These results are then used to define the response of optimized thermal conductivities to possible errors in temperature data.
In general, temperature deviations resulting from perturbation δk in the effective thermal conductivity vector k* can for time tn and depth zi be described by the equation
where represents the local temperature responses at tn and zi to unit perturbations in kj, . The full sensitivity matrix A describes the spatio-temporal distribution of the deviations between temperature simulations carried out using perturbed and not perturbed effective thermal conductivities. A is defined as
Each row corresponds to specific indices in time, , and space, and columns correspond to different k values , and in a general case J is not fixed to 6 or 8 and can be equal to M. Then the matrix form of Eq. (9) is
To derive the J columns in the sensitivity matrix A in Eq. (10) the forward model (see Eq. 3) was run J times consecutively, perturbing the individual thermal conductivities k* j, by the value Δk. The elements of A are approximated by the finite difference formula
where the entries of δkj are zero except for the jth one, which is W (m K)−1. Note that a perturbation of k* j will change all values of k(z) between zj−1 and zj+1 due to the linear interpolation for k. The full matrix A can be used to quantify deviations in the optimal k* given an assumption regarding the possible errors in measured temperature values. Since the upper and lower boundary conditions are given by and for every k, the elements and for , in A in Eq. (10) are zero. Consequently, we do not assess the effect of errors in the temperature measurements on k* through the values used to initialize and force the forward model.
Assuming that perturbations δτ in the temperature data τ will result in deviations δk of the optimized thermal conductivity values from the original estimate k*, the optimization problem in Eq. (8) now has the solution
The weight matrix W in Eq. (15) is used to select only those values of δτ that are inside the computational domain. Columns in the sensitivity matrix Aτ correspond to the different k* j nodes and individual elements express the feedback to be expected from the respective thermal conductivity value to unit perturbations at different depths and moments in time.
It is further assumed that errors in temperature measurements δτ are independent in space and constant in time and can be expressed as
where the elements of column vector δς are normally distributed random variables , , with zero mean and standard deviations si. The elements of the MN×M matrix ℰ are zero except for ℰki=1, , , . Then the error in is . The expected values of δk can be evaluated using Eqs. (15) and (16):
and the variances are found on the diagonal of the covariance matrix
where S is a square, diagonal M×M matrix with on the diagonal. In the numerical experiments for the first four time periods when nine thermistor strings were used si is taken as the time average of in Eq. (2): and for spring 2015 when only one thermistor string was used as a vertically uniform value of .
3.3.2 Depth of temperature sensors
The process of conductive heat flow is governed by the spatial temperature gradient . Uncertainty in the empirical estimates of the gradient depends not only on the accuracy and precision of the measured temperature values, but also on the possible biases in the estimates of the depths of the sensors. Real positions of the latter may differ from the ones assumed according to the original design due to the curvature of the thermistor cable in the borehole. Here possible biases in z are converted to possible biases in T using the definition of the temperature gradient and the uncertainty in k* is quantified using Aτ in Eq. (15).
The error in vertical separation of the temperature values attributed to depths zi and zi+1 is assumed to be constant in time and is denoted by δϵi+1. Then
where , δzi is the cumulative error in depth. At z1 we have the error δz1=δϵ1. If δϵi, are normally distributed random variables independent of each other and with the mean values ϵ and variances , , then . With the lower triangular M×M matrix R in
the relation between local (δϵ) and cumulative (δz) depth errors is δz=Rδϵ.
The temperature perturbation value can be expressed using the temperature gradient and the position error δzi as
with elements of the matrix Az expressing the response of the k* value in question to unit perturbations in depths of empirical temperature values (cf. with Aτ in Eq. 15). The expected value of δk is
where eM is a vector of length M with elements 1. Since stretching of the cable in a borehole is not likely, both δϵ and 𝔼[δk] values are negative rather than positive.
In line with Eq. (18), the covariance in δk is
The variances of δk are found on the diagonal of ℂov[δk],
and show how the variances are magnified to the variances in thermal conductivity through δϵ.
In the absence of empirical estimates of deviation of the real depths of sensors with respect to the designed ones, the quantification of ϵ is based on the fact that the vertical deviation of the direction of cable with thermistors is constrained by the walls of the borehole separated by diameter d=0.05 m. The uncertainty ϵ can be expressed as
where L is the typical separation between two points on the cable where it touches the walls of the borehole. In the case in which L is assumed to be equal to 1 m, ϵ is 1.3 mm.
Both the forward model and the optimization routine rely on the empirical data on the subsurface density. To assess the sensitivity of optimized thermal conductivity values k* to errors in density data ϱ we first assess the feedback of simulated errors to deviations in ρ* and then translate these results to perturbations δk using the matrix Aτ (see Eq. 15).
The sensitivity of the temperature due to perturbations δρ in ρ is computed following the logic of Eq. (9). The sensitivity matrix B is defined as
or in matrix form
with the elements in B calculated and ordered in the same manner as in A in Eqs. (10) and (12). A perturbation δρ can be interpreted as a perturbation in the temperature in the first term of Eq. (7) and by Eq. (15) it follows that the perturbation in k* will be
The second term in Eq. (7) has no influence on the solution here since we minimize over k.
3.3.4 Hessian of the objective function
where A, B, W, and γ are defined in Eqs. (11), (27), (7), and (6) and I is the identity matrix. A small perturbation of k*, ρ* will change Fτ,ϱ by δχTHδχ. The Hessian H is positive definite with positive eigenvalues and orthogonal eigenvectors when γ>0. Thus, k*, ρ* is a local minimum of Fτ,ϱ. In the direction δχ of an eigenvector with a large eigenvalue, the optimum is well defined but in a direction along an eigenvector corresponding to a small eigenvalue the uncertainty in the optimum is larger.
4.1 Measured subsurface density and temperature
The subsurface glacier profile observed in the cores consists of snow and firn with multiple ice lenses (Fig. 1). The measured snow and firn density (ϱ) varies from 350 to 900 kg m−3 with a gradual increase downward and occasional spikes corresponding to the ice layers apparent from the stratigraphical record. Compared with the measured values, the optimized densities (ρ*) show similar ranges and the overall pattern of an increase with depth (Fig. 1).
The evolution of subsurface temperature measured during five periods used for the simulations is shown in Fig. 2. The position of the upper boundary shifts upward following the snow accumulation at the surface. The evolution of deep temperature outside of the melt season generally follows the surface temperature with significant dampening with depth and time delay of the amplitude. During autumn 2014 temperature continuously decreases at all depths. During the four spring seasons the same pattern is observed only below ca. 3–4 m, while the upper part of the profile experienced warming. The measured temperature generally increases with depth; however, the opposite tendency is observed for the upper ca. 1–2 m of the profile in the spring seasons, particularly towards the end of the simulation periods. The simulation domains where Eq. (3) is solved and k and ρ are optimized in Eq. (8) are bound by the white curves limiting the area with temperature values colder than −2 ∘C and depths at the time of instrument installation larger than 1 m. The snow and firn temperature measured in April–May 2014 is significantly lower than the values registered at similar depths below the surface in April–May of the following year (Fig. 2). Based on this finding we suggest that the late part of the winter season in 2014 (February–April) was colder than in 2015. This finding is supported by the data from an automated weather station (AWS) at Nordenskiöldbreen (600 m a.s.l.), according to which in 2015 March and April were warmer than in 2014 by 2.5 and 6 ∘C correspondingly (description of AWS and data can be found at: http://www.projects.science.uu.nl/iceclimate/aws/, last access: 24 June 2019). The vertical temperature gradients measured in 2015 on 11 and 15 April by the thermistor strings installed a year earlier and by the new installation are in good correspondence (Fig. 2). This suggests that during 1 year (April 2014–April 2015) gravitational densification of snow and firn did not result in significant change of the separation between sensors and justifies the usage of time-constant density profiles in our simulations.
4.2 Optimized thermal conductivity values k*
The effective thermal conductivities k* optimized according to the Eq. (8) range from 0.2 to 1.6 W (m K)−1 and are presented in Fig. 3b. The values consistently increase with depth at the rate of ca. 0.11 W (m2 K)−1 on average and somewhat slower in the autumn of 2014. The temporal change in k* values can also be assessed with reference to Fig. 3b since the profiles for different domains are plotted with a vertical offset accounting for the surface accumulation. The overall tendency is increase in k* over time with an average rate of 0.09 W (m K yr)−1. Provided that the surface accumulation rates at the field site are slightly above 1 m, this is less than the expected value from the vertical gradient in k*. In the absence of a physical process that could result in a decrease in the firn effective thermal conductivity k* over time, the apparent lowering of optimized k* values, such as seen between 8 and 10 m from spring 2012 to spring 2013, is attributed to the uncertainties in empirical data such as temperature readings, depths of individual sensors, density measurements and estimates of the surface accumulation rate.
The uncertainty bar for each k* value is calculated following Eq. (18) and denotes the intervals of 1 standard deviation from the mean value. The bars take into account only the possible errors in subsurface temperature measurements quantified through the time-averaged standard deviations in values from different thermistor strings (Fig. 3a). Compared with the k*, the uncertainty values are generally small and are between ca. 0.03 and 0.25 W (m K)−1. The overall increase in uncertainty of k* values with depth is due to the downward decrease in temporal and spatial temperature gradients and also in the number of measurements with temperature below −2 ∘C. Altogether this results in an increase in chances for the optimization routine to converge to a “wrong” k* value based on biased temperature simulations resulting from it, as the biases are smaller and lesser in number than in the upper part of the profile. The deepest k* values in spring 2012 and 2013 exhibit outstandingly high uncertainty values, which is explained by the fact that the values are found at depths where temperature never reaches below −2 ∘C and are thus constrained only through the part of the linear fit to the k* value above that lies within the simulation domain. Due to the lack of data and large variances of the temperature measurements deep down, the deepest k* values in each of the two domains in 2012 and 2013 are likely not reliable.
In Fig. 4, the optimal k* (a) and ρ* (b) values are computed for the spring 2012 domain for different γ values. The result varies for the extreme values of γ but is consistent for intermediate γ. The k* and ρ* values are indistinguishable for γ=10, 103 and γ=10, 103, 107, respectively. The eigenvector for the smallest eigenvalue (Fig. 4c) shows that the uncertainty is largest for the deepest k* values when γ≥10 in agreement with Fig. 3b. When then the regularization of ρ is insufficient with oscillations in the solution and an eigenvector with large entries and uncertainty for all ρ values. The results for other seasons are similar.
4.3 Sensitivity of k* to errors in temperature measurements
The results of sensitivity experiments exploring the feedback (Aτ in Eq. 15) of k* values to errors in subsurface temperature measurements are presented in Fig. 5. The black markers close to the right border of each dataset show the positions of k* points. The colour at any specific point in depth zi and time tn corresponds to the response of the k* j value at the depth indicated by the black circle marker to a unit change in temperature (1 ∘C) at depth zi and time tn. The sensitivity is set to zero outside of the computational domain, the lower boundary of which is shown by the black curves.
It can also be noted that here we analyse the k* feedback to errors in individual temperature values within the dataset used to constrain the optimization routine. For the first four calculation domains (spring 2012–2014 and autumn 2014), these data are the result of lateral averaging of data from nine (q=9 in Eq. 1) thermistor strings. Although the strings are coupled to the same data logger, the errors in temperature measurements can be assumed to be at least partly independent. Thus sensitivities Aτ to temperature errors coming from single thermistor strings and not laterally averaged datasets can be expected to be lower than indicated by colour coding in Fig. 5 by a factor of .
Optimized k* values are not very sensitive to single errors in subsurface temperature data: for j=1…7 (Fig. 5a–g) the expected response of the thermal conductivity values k* j to a 1 ∘C error varies between and W (m K)−1, corresponding to the and , where and are the mean and standard deviation of the Aτj values. However, as it was demonstrated earlier (see the error bars in Fig. 3b), a systematic time-independent bias in temperature data can result in notable deviations of the k* estimates.
Sensitivities Aτ for the autumn 2014 domain exhibit a distinctively different range and pattern of spatio-temporal change of values with respect to the spring seasons. The reasons for that are not completely understood and elucidation may require additional empirical data from other autumn seasons. It can be noted that during the period from September 2014 to April 2015 the dominant tendency in the change of surface temperature was decrease, which induced continuous cooling of the subsurface profile.
Several patterns in the spatio-temporal distribution of the sensitivities Aτ can be noted. In most cases the sign of the sensitivity is positive between the depths of and k* j and negative between the depths of k* j and if k* j is the k value being tested for sensitivity (depth is marked by the black circle in Fig. 5). This pattern is reversed when considering depth levels further away from the circles (between and and between and ) and the period of alternation roughly corresponds to the spacing between k* values, which is 1 m. This result is explained through changes in the vertical temperature gradient induced by the temperature perturbations. The overall pattern in the vertical change of temperature is to increase downwards; thus a temperature increase at a certain depth results in increase in the temperature gradient just above and decrease just below that depth. These changes in temperature gradient are respectively compensated for by negative and positive deviations in k values. Due to the piecewise linear interpolation of the k profile based on J k* values, a perturbation in also affects the thermal conductivities below and above it. Therefore and adjacent to k* will display the tendency opposite to the one demonstrated by k* j to compensate for the associated changes in the heat flux. This pattern is also apparent when comparing different panels in Fig. 5.
In most cases k* is most sensitive to temperature errors ca. 0.5–1 m above and below its evaluation level, which is evidenced by the more intense colours in the sensitivity fields found in the vicinity of the black circle markers. The amplitude of cycles with alternating sign demonstrated by the sensitivity Aτ fades away with distance from the perturbed k*. The k* values found deeper down in the vertical profile appear to be more sensitive to errors in temperature, as is seen in more intensive colours around the circular markers in panels (e)–(h) of Fig. 5 compared to panels (a)–(d) in the same figure.
4.4 Sensitivity of k* to errors in depths of temperature measurements
The feedback of optimized thermal conductivities to errors in depths of temperature values used to constrain the optimization routine (Az in Eq. 21) is presented in Fig. 6, where different panels show results for the five simulation domains. The colour of a point corresponding to index j and depth z indicates the expected bias in k* j found at depth highlighted by the black circle in column j given a negative bias of −1 cm in depth of temperature values at depth z. The sensitivities Az for the lowermost k* j values (j=8 for the spring domains and j=6 for the autumn 2014) are significantly larger than for other k* j nodes due to propagation of the corresponding high Aτ values (Fig. 5g, h) to the Az matrices as is shown in Eqs. (15) and (21). These results are not shown.
In Fig. 6 blue anomalies in the vicinity of the markers indicate that the largest response of k* j can be expected to depth errors just below and just above the depth zj. This pattern is expected since the assumed depth perturbations are negative, which increases the vertical temperature gradient and is compensated for by lower k* values to keep the same heat flux. Secondly, the alternating pattern in the sensitivity of k*, noted in the previous section, can also be seen in the Az matrices, particularly for the deeper thermal conductivities with larger indices j: farther away in the vertical direction from the circular markers the negative anomalies in Az are replaced by less negative and even positive values and then switch back to the significantly negative range. This behaviour of the sensitivity is caused by the piecewise linear interpolation applied to derive the 0.1 m spaced k profile used in the forward model from 1 m spaced optimized k* values. In an attempt to preserve the heat flux and minimize the misfit between measured and simulated measurements the optimization routine will tend to overestimate the and values in case k* j is forced to be too low due to the perturbation in depths of temperature values. The uncertainties Az in Fig. 6 are notably larger above the circular markers. The thermal conductivities are much more sensitive to errors in depths of temperature values occurring above the depth level than below because of the accumulation of the position errors from upper levels to the bottom. Another reason is probably that the vertical temperature gradients are larger in the upper part of the profile. At the same time, due to the influence of the cable weight, thermistor strings can be expected to experience less coiling in the shallow part of a borehole, possibly compensating for the larger Az values there.
Assuming that thermistor strings are in contact with the borehole walls every L=1 m (see Eq. 25), the mean values of δk in Eq. (22) are on the order of −0.02 to 0.07 W (m K)−1 for k* 1–k* 6 and significantly larger lower down in the profile, where temperature gradients are lower and the number of available empirical data is less (Table 1). The expected uncertainty is 3–8 times larger for the domain covering autumn 2014, which is most likely caused by the larger temperature gradients during the early period of subsurface cooling. It can be noted that the uncertainty δk in Eq. (22) is scaled by ϵ in Eq. (25), suggesting that if L=0.5 m, the values in Table 1 will increase almost by a factor of 2. The results for the variance in Table 1 using Eq. (24) show how the variance in the positional error is magnified as a variance of the error in k. If σz≈ϵ then for most k values, a rather small standard deviation.
4.5 Sensitivity of k* to errors in density
The results of experiments exploring the sensitivity of optimized thermal conductivities to possible errors in density are presented in Fig. 7. All panels in the figure are dominated by colours corresponding to the positive sensitivity values, suggesting that the general tendency is an increase in k* in response to an overestimation of density in the empirical dataset. Negative values correspond to k* j values found in the lower part of the profile (j=6…8), which are generally less reliable.
Thermal conductivity appears to be more sensitive to density errors occurring above the depth level in question. This particularly applies to the biases in ρ2. It is still in the upper part of the subsurface profile were density is relatively low and the relative importance of a 1 kg m−3 mistake in density is large. These values are not exceeded by the Aρ1 found above, most probably because of the constraints imposed by the upper boundary condition on temperature errors: Aτ is always zero at the upper boundary of the calculation domain.
In the absence of empirical data to quantify possible errors in density, it can be noted that Aρ values presented in Fig. 7 can be scaled by any assumed value 𝔼[δρ] to derive errors 𝔼[δk] to be expected from k* values. It follows that a bias of 50 kg m−3 will result in a k* deviation of up to 0.1 W (m K)−1.
4.6 Comparison of k* with earlier published results
The relation between optimized effective thermal conductivities and densities is shown in Fig. 8 using markers, and the associated linear fit,
is illustrated by the thick black line. Also shown are the predictions from similar published fits. The Sturm et al. (1997) approximation relies on multiple needle probe measurements. The Calonne et al. (2011) and Riche and Schneebeli (2013) parameterizations are constrained by k values quantified using numerical modelling of the effective thermal conductivity tensor based on detailed three-dimensional micro-tomographic models of snow samples. It can also be noted that all three parameterizations are based on k measured in seasonal snow and the datasets included only a few samples with density above 500 kg m−3, while our results are based on measurements in firn with generally higher density.
Almost all k* values are larger than the effective thermal conductivities predicted by the first two parameterizations and the difference increases for larger densities. At the same time, the linear fit to k* and ρ* pairs appears to be closest to the Riche and Schneebeli (2013) parameterization. The latter is based on measurements made on faceted and depth hoar samples developing under a strong temperature gradient and resulting in anisotropy of the bulk thermal conductivity with larger k in the vertical direction. At Lomonosovfonna, the faceted crystals developing close to the surface are likely affected by the temperate conditions during the melt season. Below the depth of ca. 1–2 m the vertical temperature gradients are not as high and the mobility of water vapour in pores is reduced due to higher density, which altogether should limit the development of microstructural anisotropy. At the same time, it can be hypothesized that the preferential water flow in snow/firn, reported from the site (Marchenko et al., 2017b), can result in prominent vertically elongated structures below the surface of Lomonosovfonna, favouring anisotropy at a larger spatial scale and faster heat transfer in the vertical direction. The obvious variability in dependence between k* and ρ* values across different domains is a further indication of the fact that the vertical dynamics in k* are caused not only by the changes in density, and proxies describing the snow/firn structure at the scale of processes active in its metamorphism are to be included in the k=f(ρ) functions along with the density (Löwe et al., 2013). At the scale of grains such data can be derived using X-ray tomography (e.g. Kaempfer et al., 2005) and at the scale of 0.1–1 m detailed radar surveys can be used (e.g. Dunse et al., 2008; Marchenko et al., 2017a).
The findings regarding effective thermal conductivity values presented above can also be compared with the results of Giese and Hawley (2015), who applied Fourier analysis to continuous temperature measurements and derived the thermal diffusivity value (κ) of 25±3 m2 yr−1 in the top 30 m of firn pack at Summit Station, Greenland. Based on our optimized effective thermal conductivity (k*) and density (ρ*) values and the specific heat capacity of ice (C) at −10 ∘C (Cuffey and Paterson, 2010, p. 400), the thermal diffusivity (κ*) is calculated as . The resulting κ* values lie in the range from 7.55 to 48.77 m2 yr−1 with the mean value and standard deviations of 25.48 and 7.52 m2 yr−1, respectively, providing a prominent match with the results from Giese and Hawley (2015). It can be noted that quantification of k* using an optimization technique requires a much less extensive dataset in terms of time and depth coverage. Furthermore, since infiltration of meltwater in the summer interrupts the conductive heat exchange in firn at Lomonosovfonna, it is not possible to apply Fourier analysis used by Giese and Hawley (2015) on our data.
The evolution of subsurface temperature was measured in firn at Lomonosovfonna, Svalbard, using several thermistor strings during April 2012 and July 2015. The data cover five periods when the subsurface profile is at least partly at subfreezing conditions. Combined with the density measurements from four cores it was used to reconstruct the effective thermal conductivity and the density of the firn layers. For that we applied an optimization routine minimizing the mean squared difference between the measured and simulated temperature evolutions and the measured and computed density.
The optimized effective thermal conductivity k* of the firn pack at Lomonosovfonna varies from 0.4 to 1.05 W (m K)−1 increasing downwards in a maximal likelihood approach for all the time periods. According to the results of sensitivity analysis, k* is not very sensitive to systematic temperature offsets. Overestimation of the separation between sensors resulting from possible tortuosity of the cable in the borehole leads to overestimation of the k* values. Positive deviation in density estimates also results in overestimation of the k* values, while negative density biases lead to an underestimation of effective thermal conductivity.
The k* results are notably higher than the k values predicted by widely used empirical parameterizations based on the firn density measurements (Sturm et al., 1997; Calonne et al., 2011) and originally constrained by measurements of ρ and k in seasonal snow. This suggests a possible underestimation of the subsurface heat fluxes by firn models relying on the equations. In regions with a climate similar to that observed in Svalbard, this is of particular importance for the cold season as the period of conductive cooling is significantly longer than conductive warming occurring in spring before the onset of melt. Thus the refreezing capacity of the firn pack at Lomonosovfonna by the onset of melt is likely to be underestimated when simulated using the parameterizations by Sturm et al. (1997) and Calonne et al. (2011).
SM designed the optimization and experiments, collected the empirical data, and performed pilot calculations. GC developed the numerical methods and performed the numerical experiments. PL designed the uncertainty experiments and supervised the software development. VP, RP, WvP and CR acquired funding and participated in field operations. All authors participated in the writing and editing of the original draft.
The authors declare that they have no conflict of interest.
This publication is contribution 92 of the Nordic Centre of Excellence SVALI funded by the Nordic Top-level Research Initiative. Authors appreciate the constructive feedback provided by the editor and two reviewers (Henning Löwe and Edwin Waddington), their efforts helped to significantly improve the paper. Authors also acknowledge additional funding from the Ymer-80 foundation, Swedish Polar Research Secretariat, Polar Program of the Netherlands Organization for Scientific Research (NWO), Arctic Field Grant of the Research Council of Norway and the Margit Althins stipend of the Royal Swedish Academy of Sciences. Logistical support during field campaigns was provided by the Norwegian Polar Institute and the University Centre in Svalbard (UNIS).
This research has been supported by Vetenskapsrådet (grant no. 621-2014-3735) and Formas (grant no. 214-2013-1600).
This paper was edited by Martin Schneebeli and reviewed by Henning Loewe and Edwin Waddington.
Bartelt, P. and Lehning, M.: A physical SNOWPACK model for the Swiss avalanche warning: Part I: numerical model, Cold Reg. Sci. Tech., 35, 123–145, https://doi.org/10.1016/S0165-232X(02)00074-5, 2002. a
Calonne, N., Flin, F., Morin, S., Lesaffre, B., du Roscoat, S. R., and Geindreau, C.: Numerical and experimental investigations of the effective thermal conductivity of snow, Geophys. Res. Lett., 38, L23501, https://doi.org/10.1029/2011GL049234, 2011. a, b, c, d, e, f
Dunse, T., Eisen, O., Helm, V., Rack, W., Steinhage, D., and Parry, V.: Characteristics and small-scale variability of GPR signals and their relation to snow accumulation in Greenland's percolation zone, J. Glaciol., 54, 333–342, https://doi.org/10.3189/002214308784886207, 2008. a
Giese, A. L. and Hawley, R. L.: Reconstructing thermal properties of firn at Summit, Greenland, from a temperature profile time series, J. Glaciol., 61, 503–510, https://doi.org/10.3189/2015JoG14J204, 2015. a, b, c, d, e
Jordan, R., Albert, M., and Brun, E.: Physical processes within the snow cover and their parameterization, in: Snow and climate: physical processes, surface energy exchange and modeling, Cambridge University Press, Cambridge, England, 12–69, 2008. a
Ligtenberg, S. R. M., Helsen, M. M., and van den Broeke, M. R.: An improved semi-empirical model for the densification of Antarctic firn, The Cryosphere, 5, 809–819, https://doi.org/10.5194/tc-5-809-2011, 2011. a
Löwe, H., Riche, F., and Schneebeli, M.: A general treatment of snow microstructure exemplified by an improved relation for thermal conductivity, The Cryosphere, 7, 1473–1480, https://doi.org/10.5194/tc-7-1473-2013, 2013. a
Marchenko, S., Pohjola, V. A., Pettersson, R., van Pelt, W. J. J., Vega, C. P., Machguth, H., Bøggild, C. E., and Isaksson, E.: A plot-scale study of firn stratigraphy at Lomonosovfonna, Svalbard, using ice cores, borehole video and GPR surveys in 2012–14, J. Glaciol., 63, 67–78, https://doi.org/10.1017/jog.2016.118, 2017a. a, b, c, d, e
Marchenko, S., van Pelt, W. J. J., Claremar, B., Pohjola, V., Pettersson, R., Machguth, H., and Reijmer, C.: Parameterizing Deep Water Percolation Improves Subsurface Temperature Simulations by a Multilayer Firn Model, Frontiers in Earth Science, 5, 16, https://doi.org/10.3389/feart.2017.00016, 2017b. a, b, c, d
Marchenko, S., Pohjola, V., Pettersson, R., van Pelt, W., Vega, C., and Isaksson, E.: Density and stratigraphy of firn at Lomonosovfonna derived from shallow cores in 1997–2015, PANGAEA, https://doi.org/10.1594/PANGAEA.902221, 2019b. a
Morin, S., Domine, F., Arnaud, L., and Picard, G.: In-situ monitoring of the time evolution of the effective thermal conductivity of snow, Cold Reg. Sci. Technol., 64, 73–80, https://doi.org/10.1016/j.coldregions.2010.02.008, 2010. a
Nicolsky, D. J., Romanovsky, V. E., and Tipenko, G. S.: Using in-situ temperature measurements to estimate saturated soil thermal properties by solving a sequence of optimization problems, The Cryosphere, 1, 41–58, https://doi.org/10.5194/tc-1-41-2007, 2007. a
Oldroyd, H., Higgins, C., Huwald, H., Selker, J., and Parlange, M.: Thermal diffusivity of seasonal snow determined from temperature profiles, Adv. Water Resour., 55, 121–130, https://doi.org/10.1016/j.advwatres.2012.06.011, 2013. a
Osokin, N. and Sosnovsky, A.: Field investigation of efficient thermal conductivity of snow cover on Spitsbergen, Ice and Snow, 54, 50–58, https://doi.org/10.15356/2076-6734-2014-3-50-58, 2014 (in Russian). a
Pälli, A., Kohler, J. C., Isaksson, E., Moore, J. C., Pinglot, J. F., Pohjola, V. A., and Samuelsson, H.: Spatial and temporal variability of snow accumulation using ground-penetrating radar and ice cores on a Svalbard glacier, J. Glaciol., 48, 417–424, https://doi.org/10.3189/172756502781831205, 2002. a
Riche, F. and Schneebeli, M.: Thermal conductivity of snow measured by three independent methods and anisotropy considerations, The Cryosphere, 7, 217–227, https://doi.org/10.5194/tc-7-217-2013, 2013. a, b, c, d, e, f
Sergienko, O. V., MacAyeal, D. R., and Thom, J. E.: Reconstruction of snow/firn thermal diffusivities from observed temperature variation: application to iceberg C16, Ross Sea, Antarctica, 2004–07, Ann. Glaciol., 49, 91–95, https://doi.org/10.3189/172756408787814906, 2008. a, b, c
van Pelt, W. J. J., Oerlemans, J., Reijmer, C. H., Pohjola, V. A., Pettersson, R., and van Angelen, J. H.: Simulating melt, runoff and refreezing on Nordenskiöldbreen, Svalbard, using a coupled snow and energy balance model, The Cryosphere, 6, 641–659, https://doi.org/10.5194/tc-6-641-2012, 2012. a, b
van Pelt, W. J. J., Oerlemans, J., Reijmer, C. H., Pettersson, R., Pohjola, V. A., Isaksson, E., and Divine, D.: An iterative inverse method to estimate basal topography and initialize ice flow models, The Cryosphere, 7, 987–1006, https://doi.org/10.5194/tc-7-987-2013, 2013. a
van Pelt, W. J. J., Pettersson, R., Pohjola, V. A., Marchenko, S., Claremar, B., and Oerlemans, J.: Inverse estimation of snow accumulation along a radar transect on Nordenskiöldbreen, Svalbard, J. Geophys. Res.-Earth, 119, 816–835, https://doi.org/10.1002/2013JF003040, 2014. a
Weertman, J.: Creep Deformation of Ice, Annu. Rev. Earth Pl. Sc., 11, 215–240, https://doi.org/10.1146/annurev.ea.11.050183.001243, 1983. a
Wendl, I. A.: High resolution records of black carbon and other aerosol constituents from the Lomonosovfonna 2009 ice core, PhD thesis, Departement für Chemie und Biochemie der Universität Bern, 2014. a
Yang, C.-Y.: A linear inverse model for the temperature-dependent thermal conductivity determination in one-dimensional problems, Appl. Math. Model., 22, 1–9, https://doi.org/10.1016/S0307-904X(97)00101-7, 1998. a
Zhang, T. and Osterkamp, T.: Considerations in determining thermal diffusivity from temperature time series using finite difference methods, Cold Reg. Sci. Tech., 23, 333–341, https://doi.org/10.1016/0165-232X(94)00021-O, 1995. a