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)⁻¹ 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.

1 Introduction

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): $Q=-k\mathrm{\nabla }T$. Since most temperature fluctuations occur at the surface, the dominant direction of heat flux in near-surface glacier layers is vertical: $Q=-k\frac{\mathrm{d}T}{dz}$, 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 $\mathit{\kappa }=k(\mathit{\rho }C{)}^{-\mathrm{1}}$.

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)⁻¹ 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.

2 Field data

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⁻¹. 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⁻¹. 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.

Figure 1 Subsurface stratigraphy and density in four cores recovered in April 2012 (a), 2013 (b), 2014 (c) and 2015 (d). Also shown are the density profiles simulated using the maximum likelihood approach in Eq. (8).

2.2 Subsurface temperature

2.2.1 Equipment

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.

Figure 2 Measured evolution of subsurface temperature. The vertical axis is referenced to the glacier surface in April 2015. The horizontal axis is referenced to 21 April 2012 and is linear for the periods covered by measurements: 21 April–19 June 2012 (“spring 2012”), 22 April–1 June 2013 (“spring 2013”), 18 April–4 July 2014 (“spring 2014”), 25 September 2014–11 April 2015 (“autumn 2014”), and 15 April–9 July 2015 (“spring 2015”) and arbitrary between the domains. White curves indicate the upper and lower boundaries of the simulation domains. The colour bars along the horizontal axis show the intervals between measurements: yellow – 1 h, blue – 3 h and green – 12 h.

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):

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_i^n={\displaystyle \frac{\mathrm{1}}{q}}\sum _{m=\mathrm{1}}^q{\mathit{\tau }}_{im}^n,\end{array}$$

where ${\mathit{\tau }}_{im}^n$ 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

$$\begin{array}{}\text{(2)}& {\left({\mathit{\sigma }}_{\mathit{\tau }i}^n\right)}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{q-\mathrm{1}}}\sum _{m=\mathrm{1}}^q{\left({\mathit{\tau }}_{im}^n-{\mathit{\tau }}_i^n\right)}^{\mathrm{2}}.\end{array}$$

3 Modelling

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

$$\begin{array}{}\text{(3)}& \mathit{\rho }\left(z\right)C{\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(k\right(z\left){\displaystyle \frac{\partial T}{\partial z}}\right),\end{array}$$

where ρ is the subsurface density (kg m⁻³), C=2027 J (kg K)⁻¹ 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)⁻¹). 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 $T_i^n$ be the temperature T(z_i, tⁿ) at tⁿ=nΔt, $n=\mathrm{1},\mathrm{\dots },N$ and z_i=iΔz, $i=\mathrm{1},\mathrm{\dots },M^n$. The number of time steps is N and the number of nodes in space Mⁿ varies in time. The solution at z_i is advanced in time from tⁿ to tⁿ⁺¹ as

$$\begin{array}{}\text{(4)}& {\mathit{\rho }}_{i}C{\displaystyle \frac{T_i^{n+\mathrm{1}}-T_i^n}{\mathrm{\Delta }t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathbf{D}{T}_i^{n+\mathrm{1}}+\mathbf{D}{T}_i^n\right),n=\mathrm{1},\mathrm{\dots },N-\mathrm{1},\end{array}$$

where $\mathbf{D}{T}_i^n$ is the spatial temperature derivative at z_i approximated by

$$\begin{array}{}\text{(5)}& \mathbf{D}{T}_i^n={\displaystyle \frac{\left(k_{i+\mathrm{1}}+k_i\right)\left(T_{i+\mathrm{1}}^n-T_i^n\right)-\left(k_i+k_{i-\mathrm{1}}\right)\left(T_i^n-T_{i-\mathrm{1}}^n\right)}{\mathrm{2}\mathrm{\Delta }{z}^{\mathrm{2}}}}.\end{array}$$

Collecting terms in Eq. (4) with Tⁿ⁺¹ on the left-hand side and terms with Tⁿ on the right-hand side, a tridiagonal system of linear equations can de derived and solved for $T_i^{n+\mathrm{1}}$ 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 Mⁿ nodes from the top at z₁=1 m to the bottom at $z_{M^n}$ 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 $z_{M^n}$ 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 ${\mathit{\tau }}_i^{\mathrm{1}}$. The upper and lower boundary conditions are of Dirichlet type and are determined by the temperatures measured at z₁ and $z_{M^n}$: $T_{\mathrm{1}}^n={\mathit{\tau }}_{\mathrm{1}}^n$ and $T_{M^n}^n={\mathit{\tau }}_{M^n}^n$.

The density ρ(z) and effective thermal conductivity k(z) at depths z_i 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 $F_{\mathit{\tau },\mathit{\varrho }}(k,\mathit{\rho })$. Following Smith (2013), the latter is defined by the sums of squared differences between the simulated and measured temperature and density:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{F}_{\mathit{\tau },\mathit{\varrho }}(k,\mathit{\rho })& =\sum _{n=\mathrm{1}}^N\sum _{i=\mathrm{1}}^{M^n}{\displaystyle \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{{\mathit{\tau }}_i^n}\right)}^{\mathrm{2}}}}{\left(T_i^n(k,\mathit{\rho })-{\mathit{\tau }}_i^n\right)}^{\mathrm{2}}\\ & +\mathit{\gamma }\sum _{i=\mathrm{1}}^{M^n}{\displaystyle \frac{\mathrm{1}}{{\left({\mathit{\sigma }}_{{\mathit{\varrho }}_i}\right)}^{\mathrm{2}}}}{\left({\mathit{\rho }}_i-{\mathit{\varrho }}_i\right)}^{\mathrm{2}},\end{array}\end{array}$$

and the optimization routine attempts to minimize $F_{\mathit{\tau },\mathit{\varrho }}(k,\mathit{\rho })$ by adjusting the conductivity k and density ρ values used in simulations. In Eq. (6) $T_i^n(k,\mathit{\rho })$ and ${\mathit{\tau }}_i^n$ are the simulated and measured temperature values at time tⁿ and depth z_i, and ϱ_i is the density measured at depth z_i. 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 $({\mathit{\sigma }}_{{\mathit{\tau }}_i^n}{)}^{\mathrm{2}}$ 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 $({\mathit{\sigma }}_{{\mathit{\tau }}_i^n}{)}^{\mathrm{2}}$. 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 ${\mathit{\sigma }}_{{\mathit{\varrho }}_i}$ was set to 1 kg m⁻³. 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 l²-norm. The choice of γ depends on the size of the data N, Mⁿ, 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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{T}={\left(T_{\mathrm{1}}^{\mathrm{1}},T_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },T_{\mathrm{1}}^N,T_{\mathrm{2}}^{\mathrm{1}},\mathrm{\dots },T_{\mathrm{2}}^N,T_{\mathrm{3}}^{\mathrm{1}},\mathrm{\dots },T_M^N\right)}^T,\\ {\displaystyle }& {\displaystyle }\mathit{\tau }={\left({\mathit{\tau }}_{\mathrm{1}}^{\mathrm{1}},{\mathit{\tau }}_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },{\mathit{\tau }}_{\mathrm{1}}^N,{\mathit{\tau }}_{\mathrm{2}}^{\mathrm{1}},\mathrm{\dots },{\mathit{\tau }}_{\mathrm{2}}^N,{\mathit{\tau }}_{\mathrm{3}}^{\mathrm{1}},\mathrm{\dots },{\mathit{\tau }}_M^N\right)}^T,\\ {\displaystyle }& {\displaystyle }\mathit{\rho }={\left({\mathit{\rho }}_{\mathrm{1}},{\mathit{\rho }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\rho }}_M\right)}^T,\\ {\displaystyle }& {\displaystyle }\mathit{\varrho }={\left({\mathit{\varrho }}_{\mathrm{1}},{\mathit{\varrho }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\varrho }}_M\right)}^T,\end{array}$$

where $M={max}_n{M}^n$. Let the diagonal matrices W have $W_{jj}=\mathrm{1}/{\left({\mathit{\sigma }}_{{\mathit{\tau }}_i^n}\right)}^{\mathrm{2}}$, $j=N(i-\mathrm{1})+n$, $n=\mathrm{1},\mathrm{\dots },N$, $i=\mathrm{1},\mathrm{\dots },M^n$, and W_jj=0, $j=N(i-\mathrm{1})+n$, $n=\mathrm{1},\mathrm{\dots },N$, $i=M^n+\mathrm{1},\mathrm{\dots },M$, on the diagonal. Then F in Eq. (6) can be written as

$$\begin{array}{}\text{(7)}& F_{\mathit{\tau },\mathit{\varrho }}(k,\mathit{\rho })=(T-\mathit{\tau }{)}^T\mathbf{W}(T-\mathit{\tau })+\mathit{\gamma }(\mathit{\rho }-\mathit{\varrho }{)}^T(\mathit{\rho }-\mathit{\varrho }).\end{array}$$

The diagonal elements of W vanish when i>Mⁿ since the sums in Eq. (6) are restricted to i≤Mⁿ.

The optimal k_* and ρ_* minimize F in Eqs. (6) and (7) and the solution of the nonlinear least-squares optimization problem can be written as

$$\begin{array}{}\text{(8)}& \left({\mathit{k}}_{*},{\mathit{\rho }}_{*}\right)=\underset{\mathit{k},\mathit{\rho }}{\text{arg\hspace{0.17em}min}}{F}_{\mathit{\tau },\mathit{\varrho }}(\mathit{k},\mathit{\rho }).\end{array}$$

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.

3.3.1 Temperature

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 tⁿ and depth z_i be described by the equation

$$\begin{array}{}\text{(9)}& T_i^n\left({\mathit{k}}_{*}+\mathit{\delta }\mathit{k},{\mathit{\rho }}_{*}\right)-T_i^n\left({\mathit{k}}_{*},{\mathit{\rho }}_{*}\right)=\sum _{j=\mathrm{1}}^J{A}_{ij}^n\mathit{\delta }{k}_j,\end{array}$$

where $A_{ij}^n$ represents the local temperature responses at tⁿ and z_i to unit perturbations in k_j, $j=\mathrm{1},\mathrm{\dots },J$. 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

$$\begin{array}{}\text{(10)}& \mathbf{A}=\begin{array}{cccc}{A}_{\mathrm{11}}^{\mathrm{1}}& A_{\mathrm{12}}^{\mathrm{1}}& \mathrm{\cdots }& A_{\mathrm{1}J}^{\mathrm{1}}\\ \mathrm{⋮}& & & \mathrm{⋮}\\ A_{\mathrm{11}}^N& A_{\mathrm{12}}^N& \mathrm{\cdots }& A_{\mathrm{1}J}^N\\ & & & \\ A_{\mathrm{21}}^{\mathrm{1}}& A_{\mathrm{22}}^{\mathrm{1}}& \mathrm{\cdots }& A_{\mathrm{2}J}^{\mathrm{1}}\\ \mathrm{⋮}& & & \mathrm{⋮}\\ A_{M\mathrm{1}}^N& A_{M\mathrm{2}}^N& \mathrm{\cdots }& A_{MJ}^N\end{array}.\end{array}$$

Each row corresponds to specific indices in time, $n=\mathrm{1},\mathrm{\dots },N$, and space, $i=\mathrm{1},\mathrm{\dots },M$ and columns correspond to different k values $j=\mathrm{1},\mathrm{\dots },J$, 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

$$\begin{array}{}\text{(11)}& T\left({\mathit{k}}_{*}+\mathit{\delta }\mathit{k},{\mathit{\rho }}_{*}\right)=T\left({\mathit{k}}_{*},{\mathit{\rho }}_{*}\right)+\mathbf{A}\mathit{\delta }\mathit{k}.\end{array}$$

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, $j=\mathrm{1},\mathrm{2}\mathrm{\dots },J$ by the value Δk. The elements of A are approximated by the finite difference formula

$$\begin{array}{}\text{(12)}& A_{ij}^n={\displaystyle \frac{T_i^n\left({\mathit{k}}_{*}+\mathit{\delta }{k}_j,{\mathit{\rho }}_{*}\right)-T_i^n\left({\mathit{k}}_{*},{\mathit{\rho }}_{*}\right)}{\mathrm{\Delta }k}},\end{array}$$

where the entries of δk_j are zero except for the jth one, which is $\mathrm{\Delta }k={\mathrm{10}}^{-\mathrm{6}}$ W (m K)⁻¹. Note that a perturbation of k_* j will change all values of k(z) between z_j−1 and z_j+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 ${\mathit{\tau }}_{\mathrm{1}}^n$ and ${\mathit{\tau }}_{M^n}^n$ for every k, the elements $A_{\mathrm{1}j}^n$ and $A_{M^{n}j}^n$ for $j=\mathrm{1},\mathrm{\dots },J$, $n=\mathrm{1},\mathrm{\dots },N$ 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

$$\begin{array}{}\text{(13)}& {\mathit{k}}_{*}+\mathit{\delta }\mathit{k}={\text{arg\hspace{0.17em}min}}_k{F}_{\mathit{\tau }+\mathit{\delta }\mathit{\tau },\mathit{\varrho }}\left(\mathit{k},{\mathit{\rho }}_{*}\right).\end{array}$$

If δτ in Eq. (13) is small and ρ_* is fixed, then using Eqs. (7) and (11) one can show that δk in Eq. (13) satisfies

$$\begin{array}{}\text{(14)}& {\mathbf{A}}^T\mathbf{WA}\mathit{\delta }k={\mathbf{A}}^T\mathbf{W}\mathit{\delta }\mathit{\tau }.\end{array}$$

By linearizing the dependence of $T_j\left({\mathit{k}}_{*}+\mathit{\delta }\mathit{k},{\mathit{\rho }}_{*}\right)-T\left({\mathit{k}}_{*},{\mathit{\rho }}_{*}\right)$ on δk in Eq. (11), we arrive at a linear least-squares problem to solve for δk in Eq. (14):

$$\begin{array}{}\text{(15)}& \mathit{\delta }\mathit{k}=({\mathbf{A}}^T\mathbf{WA}{)}^{-\mathrm{1}}{\mathbf{A}}^T\mathbf{W}\mathit{\delta }\mathit{\tau }={\mathbf{A}}_{\mathit{\tau }}\mathit{\delta }\mathit{\tau }.\end{array}$$

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

$$\begin{array}{}\text{(16)}& \mathit{\delta }\mathit{\tau }=\mathcal{E}\mathit{\delta }\mathit{\varsigma },\end{array}$$

where the elements of column vector δς are normally distributed random variables $\mathit{\delta }{\mathit{\varsigma }}_i\sim \mathcal{N}(\mathrm{0},s_i^{\mathrm{2}})$, $i=\mathrm{1},\mathrm{\dots },M$, with zero mean and standard deviations s_i. The elements of the MN×M matrix ℰ are zero except for ℰ_ki=1, $k=N(i-\mathrm{1})+n$, $n=\mathrm{1},\mathrm{\dots },N$, $i=\mathrm{1},\mathrm{\dots },M$. Then the error in ${\mathit{\tau }}_i^n$ is $\mathit{\delta }{\mathit{\tau }}_i^n=\mathit{\delta }{\mathit{\varsigma }}_i$. The expected values of δk can be evaluated using Eqs. (15) and (16):

$$\begin{array}{}\text{(17)}& \mathbb{E}\left[\mathit{\delta }\mathit{k}\right]={\mathbf{A}}_{\mathit{\tau }}\mathcal{E}\mathbb{E}\left[\mathit{\delta }\mathit{\varsigma }\right]=\mathrm{0},\end{array}$$

and the variances are found on the diagonal of the covariance matrix

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}\mathbb{C}\mathrm{ov}\left[\mathit{\delta }\mathit{k}\right]& =\mathbb{E}\left[\left({\mathbf{A}}_{\mathit{\tau }}\mathit{\delta }\mathit{\tau }\right){\left({\mathbf{A}}_{\mathit{\tau }}\mathit{\delta }\mathit{\tau }\right)}^T\right]={\mathbf{A}}_{\mathit{\tau }}\mathbb{E}\left[\mathit{\delta }\mathit{\tau }\mathit{\delta }{\mathit{\tau }}^T\right]{\mathbf{A}}_{\mathit{\tau }}^T\\ & ={\mathbf{A}}_{\mathit{\tau }}\mathcal{E}\mathbb{E}\left[\mathit{\delta }\mathit{\varsigma }\mathit{\delta }{\mathit{\varsigma }}^T\right]{\mathcal{E}}^T{\mathbf{A}}_{\mathit{\tau }}^T={\mathbf{A}}_{\mathit{\tau }}\mathcal{E}\mathbf{S}{\mathcal{E}}^T{\mathbf{A}}_{\mathit{\tau }}^T,\end{array}\end{array}$$

where S is a square, diagonal M×M matrix with $s_i^{\mathrm{2}}$ on the diagonal. In the numerical experiments for the first four time periods when nine thermistor strings were used s_i is taken as the time average of ${\mathit{\sigma }}_{{\mathit{\tau }}_i^n}$ in Eq. (2): $s_i^{\mathrm{2}}=\frac{\mathrm{1}}{N}{\sum }_{n=\mathrm{1}}^N{\mathit{\sigma }}_{{\mathit{\tau }}_i^n}^{\mathrm{2}}$ and for spring 2015 when only one thermistor string was used as a vertically uniform value of $s_i^{\mathrm{2}}=\mathrm{0.02}$.

3.3.2 Depth of temperature sensors

The process of conductive heat flow is governed by the spatial temperature gradient $\frac{\partial T}{\partial z}$. 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 z_i and z_i+1 is assumed to be constant in time and is denoted by δϵ_i+1. Then

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{z}_{i+\mathrm{1}}& =z_i+\mathrm{\Delta }z+\mathit{\delta }{\mathit{ϵ}}_{i+\mathrm{1}}=z_{\mathrm{1}}+i\mathrm{\Delta }z+\sum _{j=\mathrm{1}}^{i+\mathrm{1}}\mathit{\delta }{\mathit{ϵ}}_j\\ & =z_{\mathrm{1}}+i\mathrm{\Delta }z+\mathit{\delta }{z}_{i+\mathrm{1}},\end{array}\end{array}$$

where $i=\mathrm{1},\mathrm{\dots },M^n-\mathrm{1}$, δz_i is the cumulative error in depth. At z₁ we have the error δz₁=δϵ₁. If δϵ_i, $i=\mathrm{1},\mathrm{\dots },M^n$ are normally distributed random variables independent of each other and with the mean values ϵ and variances ${\mathit{\sigma }}_z^{\mathrm{2}}$, $\mathit{\delta }{\mathit{ϵ}}_i\sim \mathcal{N}(\mathit{ϵ},{\mathit{\sigma }}_z^{\mathrm{2}})$, then $\mathit{\delta }{z}_i\sim \mathcal{N}(i\mathit{ϵ},i{\mathit{\sigma }}_z^{\mathrm{2}})$. With the lower triangular M×M matrix R in

$${\displaystyle }\mathbf{R}=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{1}& \mathrm{1}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}\\ \mathrm{1}& \mathrm{1}& \mathrm{\cdots }& \mathrm{1}\end{array}\right],$$

the relation between local (δϵ) and cumulative (δz) depth errors is δz=Rδϵ.

The temperature perturbation value $\mathit{\delta }{\mathit{\tau }}_i^n$ can be expressed using the temperature gradient and the position error δz_i as

$$\begin{array}{}\text{(20)}& \mathit{\delta }{\mathit{\tau }}_i^n=\mathit{\delta }{z}_i{\left.{\displaystyle \frac{\partial T}{\partial z}}\right|}_{z=z_i,t=t^n}.\end{array}$$

Let the elements of the matrix 𝒯 be zero except for ${\mathcal{T}}_{ki}={\left.\frac{\partial T}{\partial z}\right|}_{z=z_i,t=t^n}$, where $k=N(i-\mathrm{1})+n$, $n=\mathrm{1},\mathrm{\dots },N$, $i=\mathrm{1},\mathrm{\dots }M$. Then $\mathit{\delta }\mathit{\tau }=\mathcal{T}\mathit{\delta }z=\mathcal{T}\mathbf{R}\mathit{\delta }\mathit{ϵ}$. Thus, the minimal δk in Eq. (13) satisfies Eq. (14) with the solution

$$\begin{array}{}\text{(21)}& \mathit{\delta }\mathit{k}=({\mathbf{A}}^T\mathbf{WA}{)}^{-\mathrm{1}}{\mathbf{A}}^T\mathbf{W}\mathcal{T}\mathbf{R}\mathit{\delta }\mathit{ϵ}={\mathbf{A}}_z\mathit{\delta }\mathit{ϵ},\end{array}$$

with elements of the matrix A_z 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

$$\begin{array}{}\text{(22)}& \mathbb{E}\left[\mathit{\delta }\mathit{k}\right]={\mathbf{A}}_z\mathbb{E}\left[\mathit{\delta }\mathit{ϵ}\right]=\mathit{ϵ}{\mathbf{A}}_z{\mathit{e}}_M,\end{array}$$

where e_M 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

$$\begin{array}{}\text{(23)}& \mathbb{C}\mathrm{ov}\left[\mathit{\delta }\mathit{k}\right]={\mathbf{A}}_z\mathbb{C}\mathrm{ov}\left[\mathit{\delta }\mathit{ϵ}\right]{\mathbf{A}}_z^T={\mathit{\sigma }}_z^{\mathrm{2}}{\mathbf{A}}_z{\mathbf{A}}_z^T.\end{array}$$

The variances of δk are found on the diagonal of ℂov[δk],

$$\begin{array}{}\text{(24)}& \mathbb{V}\mathrm{ar}\left[\mathit{\delta }{k}_j\right]={\mathit{\sigma }}_z^{\mathrm{2}}\sum _{i=\mathrm{1}}^{M^n}{\mathbf{A}}_{zji}^{\mathrm{2}},\end{array}$$

and show how the variances ${\mathit{\sigma }}_z^{\mathrm{2}}$ 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

$$\begin{array}{}\text{(25)}& \mathit{ϵ}=L-\sqrt{L^{\mathrm{2}}-d^{\mathrm{2}}},\end{array}$$

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.

3.3.3 Density

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

$$\begin{array}{}\text{(26)}& T_i^n\left({\mathit{k}}_{*},{\mathit{\rho }}_{*}+\mathit{\delta }\mathit{\rho }\right)=T_i^n\left({\mathit{k}}_{*},{\mathit{\rho }}_{*}\right)+\sum _{j=\mathrm{1}}^J{B}_{ij}^n\mathit{\delta }{\mathit{\rho }}_j,\end{array}$$

or in matrix form

$$\begin{array}{}\text{(27)}& T\left({\mathit{k}}_{*},{\mathit{\rho }}_{*}+\mathit{\delta }\mathit{\rho }\right)=T\left({\mathit{k}}_{*},{\mathit{\rho }}_{*}\right)+\mathbf{B}\mathit{\delta }\mathit{\rho },\end{array}$$

with the elements $B_{ij}^n$ 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 $\mathit{\delta }\mathit{\tau }=-\mathbf{B}\mathit{\delta }\mathit{\rho }$ in the first term of Eq. (7) and by Eq. (15) it follows that the perturbation in k_* will be

$$\begin{array}{}\text{(28)}& \mathit{\delta }k=-({\mathbf{A}}^T\mathbf{WA}{)}^{-\mathrm{1}}{\mathbf{A}}^T\mathbf{WB}\mathit{\delta }\mathit{\rho }=-{\mathbf{A}}_{\mathit{\rho }}\mathit{\delta }\mathit{\rho }.\end{array}$$

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

The behaviour of the objective function $F_{\mathit{\tau },\mathit{\varrho }}(k,\mathit{\rho })$ in Eqs. (6) and (7) close to the optimum k_*, ρ_* is determined by the Hessian matrix H:

$$\begin{array}{}\text{(29)}& \mathbf{H}=\left(\begin{array}{cc}{\mathbf{A}}^T\mathbf{WA}& {\mathbf{A}}^T\mathbf{WB}\\ {\mathbf{B}}^T\mathbf{WA}& {\mathbf{B}}^T\mathbf{WB}+\mathit{\gamma }\mathbf{I}\end{array}\right),\end{array}$$

where A, B, W, and γ are defined in Eqs. (11), (27), (7), and (6) and I is the identity matrix. A small perturbation $\mathit{\delta }\mathit{\chi }=\left(\mathit{\delta }k,\mathit{\delta }\mathit{\rho }\right)$ 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 Results and discussion

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⁻³ 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 ${\mathit{\tau }}_i^n$ 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)⁻¹ and are presented in Fig. 3b. The values consistently increase with depth at the rate of ca. 0.11 W (m² K)⁻¹ 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)⁻¹. 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.

Figure 3 (a) Time-averaged standard deviations in subsurface temperature measured by different thermistor strings (s_i). (b) Optimized effective thermal conductivity values (k_*) calculated following Eq. (8) with error bars according to ℂov[δk] in Eq. (18) assuming temperature errors shown in panel a.

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)⁻¹. 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, 10³ and γ=10, 10³, 10⁷, 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 $\mathit{\gamma }={\mathrm{10}}^{-\mathrm{5}}$ 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.

Figure 4 The optimal solutions k_* (a) and ρ_* (b) in the spring 2012 domain for different values of γ in Eq. (7). The black curve in (b) shows measured density values. (c) Absolute values of the eigenvector with the smallest eigenvalue of the Hessian H in Eq. (29).

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 ${\mathit{\tau }}_i^n$ 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 z_i and time tⁿ 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 z_i and time tⁿ. The sensitivity is set to zero outside of the computational domain, the lower boundary of which is shown by the black curves.

Figure 5 Sensitivity (A_τ in Eq. 15) of effective thermal conductivities (k_* j) to errors in the temperature measurements. Panels (a–h) correspond to $j=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{8}$, in each case the depth of k_* j is highlighted by black circle; depths of other k_* values are shown by black dots. The colour at a certain point in depth and time corresponds to the feedback of the k_* j value found at the depth shown by the black circle to a unit error in temperature (+1 ^∘C) at that depth and time. Black curves indicate the lower boundaries of the computational domains; below them the sensitivity is set to 0. The horizontal axes are the same as in Fig. 2.

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 $\sqrt{q}=\mathrm{3}$.

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 $-\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{3}}$ W (m K)⁻¹, corresponding to the ${\stackrel{\mathrm{‾}}{\mathbf{A}}}_{\mathit{\tau }j}-\mathrm{3}{\mathit{\sigma }}_{{\mathbf{A}}_{\mathit{\tau }j}}$ and ${\stackrel{\mathrm{‾}}{\mathbf{A}}}_{\mathit{\tau }j}+\mathrm{3}{\mathit{\sigma }}_{{\mathbf{A}}_{\mathit{\tau }j}}$, where ${\stackrel{\mathrm{‾}}{\mathbf{A}}}_{\mathit{\tau }j}$ and ${\mathit{\sigma }}_{{\mathbf{A}}_{\mathit{\tau }j}}$ 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 $k_{*j-\mathrm{1}}$ and k_* j and negative between the depths of k_* j and $k_{*j+\mathrm{1}}$ 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 $k_{*j-\mathrm{2}}$ and $k_{*j-\mathrm{1}}$ and between $k_{*j+\mathrm{1}}$ and $k_{*j+\mathrm{2}}$) 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 $k_{{*}_j}$ also affects the thermal conductivities below and above it. Therefore $k_{*j-\mathrm{1}}$ and $k_{*j+\mathrm{1}}$ 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 (A_z 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 A_z 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 A_z matrices as is shown in Eqs. (15) and (21). These results are not shown.

Figure 6 Sensitivity (A_z in Eq. (21)) of effective thermal conductivities k_* j to errors in depths of temperature values. Panels correspond to different calculation domains: spring 2012 (a), spring 2013 (b), spring 2014 (c), autumn 2014 (d) and spring 2015 (e). In all panels the colour of a point in column j and depth z shows the feedback on effective thermal conductivity value k_* j found at the depth highlighted by the black (white in d) circle to an error in depth of temperature value (1 cm) made at depth z.

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 z_j. 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 A_z 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 A_z 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 $k_{*j+\mathrm{1}}$ and $k_{*j-\mathrm{1}}$ values in case k_* j is forced to be too low due to the perturbation in depths of temperature values. The uncertainties A_z 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 A_z 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)⁻¹ 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 ${\mathit{\sigma }}_z^{\mathrm{2}}$ in the positional error is magnified as a variance of the error in k. If σ_z≈ϵ then $\sqrt{\mathbb{V}\mathrm{ar}\left[\mathit{\delta }k\right]}\approx \mathrm{0.01}$ for most k values, a rather small standard deviation.

Table 1 Expected deviations of thermal conductivities (𝔼[δk] in Eq. 22) and the corresponding standard deviations ($\sqrt{\mathbb{V}\mathrm{ar}\left[\mathit{\delta }k\right]}$; see Eqs. 23, 24) given the magnitude of thermistor string coiling in the boreholes is expressed by L=1 m (see Eq. 25) and σ_z=1 cm.

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.

Figure 7 Sensitivity (A_ρ in Eq. 28) of effective thermal conductivities k_* j to errors in density. Panels correspond to different calculation domains: spring 2012 (a), spring 2013 (b), spring 2014 (c), autumn 2014 (d) and spring 2015 (e). In all panels the colour of a point in column j and depth z shows the feedback of effective thermal conductivity value k_* j found at the depth highlighted by the white circle to a unit error in density (1 kg m⁻³) made at depth z.

Thermal conductivity appears to be more sensitive to density errors occurring above the depth level in question. This particularly applies to the biases in ρ₂. 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⁻³ 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⁻³ will result in a k_* deviation of up to 0.1 W (m K)⁻¹.

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,

$$\begin{array}{}\text{(30)}& k=\mathrm{0.301}\cdot {\mathrm{10}}^{-\mathrm{2}}\mathit{\rho }-\mathrm{0.724},\end{array}$$

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⁻³, while our results are based on measurements in firn with generally higher density.

Figure 8 Relation between effective thermal conductivities k and densities ρ. Markers illustrate the results of the optimization routine; the black line shows the linear fit based on all k_*–ρ_* pairs. The light blue, ruby and green curves show quadratic k=f(ρ) functions according to Sturm et al. (1997), Calonne et al. (2011), and Riche and Schneebeli (2013).

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 m² yr⁻¹ 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 ${\mathit{\kappa }}_{*}=k_{*}({\mathit{\rho }}_{*}C{)}^{-\mathrm{1}}$. The resulting κ_* values lie in the range from 7.55 to 48.77 m² yr⁻¹ with the mean value and standard deviations of 25.48 and 7.52 m² yr⁻¹, 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.

5 Conclusions

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)⁻¹ 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).