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. 5 The optimized k∗ values lie 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.

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, page 403), according to which the heat flux (Q) is proportional to its thermal conductivity (k) and to the spatial temperature gradient (∇T ): Q = −k∇T . Since most temperature fluctuations 5 occur at the surface, the dominant direction of heat flux in near surface glacier layers is vertical: Q = −k dT 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 vapor. 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 10 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 κ = k (ρ C) −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 to derive κ and k values either using Fourier-type analysis or optimization techniques. Diurnal temperature 15 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 20 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 snow pack on a drifting iceberg and Oldroyd et al. (2013) for the seasonal snow pack on an Alpine glacier.
Another option is to induce a heat flux in the snow pack 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 in snow pack itself. The former 25 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 30 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 to obtain the k tensor for relatively small snow samples. 35 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 interconnectivity and 5 anisotropy of the k property.
The purpose of the present study is to reconstruct the values of effective thermal conductivity of a thick snow 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 is has 10 not been applied for the purpose so far, although Sergienko et al. (2008) employed similar routines for a seasonal snow pack and Nicolsky et al. (2007) for permafrost. Fourier analysis applied earlier for a thick firn pack at the Summit of the Greenland ice sheet (Giese and Hawley, 2015) can not be used here due to the influence of melt water and retrieving snow 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 15 the applied models.

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 asl, which is well above the equilibrium line, estimated to be at ca 720 m asl (Van Pelt et al., 2012). The local glacier thickness is 192±5.1 m (Pettersson, 2009, unpublished dataset;20 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)

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 (Figure 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 5 2.2.1 Equipment Subsurface temperature was measured using multiple thermistors grouped in several strings. They were placed in boreholes with 5.5 cm diameter drilled by a Kovacs auger. The boreholes were then backfilled by 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 neighboring boreholes (Marchenko et al., 2017a, b). In April 2015 only one thermistor string was 10 installed.
All thermistor strings were custom manufactured at Uppsala University using a multi-leaded cable with PVC jacket and unicurve NTC thermistors. In 2015 the cable was placed in multiple 2 m long rigid plastic tubes to ensure a precise and constant separation between neighboring sensors. The sensors were fed through holes in the tubes and fixed at their outer surface for a better contact with the sounded environment. 15 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 thermistor manufacturer.

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 30 measured during July 1-September 1. 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 neighboring 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 5 τ describing subsurface temperature changes in time (superscript n) and depth (subscript i): where τ n im is the temperature measured by the m-th 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 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, page 403).
Configuring the computational domains to minimize the influence of non-conductive heat fluxes, we compute the effective 15 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 . Furthermore, temperature 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.

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 10 dependent function from (Cuffey and Paterson, 2010, page 400) for the temperature of -10 • C, and k is the effective thermal conductivity (W (m K) −1 ). Given k(z), ρ(z) along with initial and boundary conditions, Eq. (3) is solved forward in time for The numerical solution of the 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 15 differences in space and trapezoidal rule in time. Let T n i be the temperature T (z i , t n ) at t n = n∆t, n = 1, . . . , N and z i = i∆z, i = 1, . . . , M n . The number of time steps is N and the number of nodes in space M n varies in time. The solution at z i is advanced in time from t n to t n+1 as: where DT n i is the spatial temperature derivative at z i approximated by: Collecting terms in (4) with T n+1 on the left hand side and terms with T n on the right hand side, a tridiagonal system of linear equations can de derived and solved for T n+1 i 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 n nodes from the top at z 1 = 1 m to the bottom at z M n where the 25 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 Figure 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 30 using measured temperature values τ 1 i . The upper and lower boundary conditions are of Dirichlet type and are determined by the temperatures measured at z 1 and z M n : T n 1 = τ n 1 and T n M n = τ n M 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 section 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 5 and later by the results of the previous inversion iteration. For the domains covering spring seasons J equals 8 and is 6 in fall 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 amount of data to constrain a very detailed k * profile (high J). Too large J will result in an oscillatory optimal solution for k.

Inverse routine 10
The effective thermal conductivity k in (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 τ, (k, ρ). Following Smith (2013), the latter is defined by the sums of squared differences between the simulated and measured temperature and density: 15 and the optimization routine attempts to minimize F τ, (k, ρ) by adjusting the conductivity k and density ρ values used in simulations. In eq. (6) T n i (k, ρ) and τ n i are the simulated and measured temperature values at time t n and depth z i , 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 fall 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 (σ τ n i ) 2 in temperature values from 20 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 (σ τ n i ) 2 . Since only one shallow core was drilled every year, empirical data are not available to quantify the errors in density measurements, the weighting term σ i 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 25 measurements in l 2 -norm. The choice of γ depends on the size of the data N , M n as well as the magnitudes of τ and .
Choice of the cost function in (6)) does not assume any correlation between k and ρ. This relation is derived later in Section 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 ca 10 m covered by 30 thermistor measurements. The second term on the right hand side of equation (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).
Introduce the vectors: . . , M n , and W jj = 0, j = N (i − 1) + n, n = 1, . . . , N, i = M n + 1, . . . , M, on the diagonal. Then F in Eq. (6) can be written as: The diagonal elements of W vanish when i > M n since the sums in Eq. (6) are restricted to i ≤ M n .
The optimal k * and ρ * minimize F in Eq. (6) and (7) and the solution of the nonlinear least squares optimization problem can be written as: It is solved by the lsqnonlin function in MATLAB.

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 in (Smith, 2013). The results of the sensitivity estimates are only valid for relatively small errors on the order of ca 15 5 -10 % of the parameter value in question.

Temperature
Errors δτ in the measured temperature values τ in Eq. (6)  In general, temperature deviations resulting from perturbation δk in the effective thermal conductivity vector k * can for time t n and depth z i be described by the equation: where A n ij represent the local temperature responses at t n and z i to unit perturbations in k j , j = 1, . . . , J. The full sensitivity 25 matrix A describes the spatio-temporal distribution of the deviations between temperature simulations done using perturbed and not perturbed effective thermal conductivities. A is defined as: Each row corresponds to specific indices in time, n = 1, . . . , N , and space, i = 1, . . . , M and columns correspond to different k values j = 1, . . . , 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: 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 = 1, 2 . . . , J by the value ∆k. The elements of A are approximated by the finite difference formula: where the entries of δk j are zero except for the j:th one which is ∆k = 10 −6 W (m K) −1 . 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 τ n 1 and τ n M n for every k, the elements A n 1j and A n M n j for j = 1, . . . , J, n = 1, . . . , 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. 15 If we now assume that δτ are perturbations in the temperature data τ , that 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: If δτ in Eq. (13) is small and ρ * is fixed, then using Eq. (7) and Eq. (11) one can show that δk in Eq. (13) satisfies: 20 By linearising the dependence of T j (k * + δk, ρ * ) − T (k * , ρ * ) on δk in Eq. (11), we arrive at a linear least squares problem to solve for δk in Eq. (14): 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 element 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: and the variances are found on the diagonal of the covariance matrix

Depth of temperature sensors
The process of conductive heat flow is governed by the spatial temperature gradient ∂T ∂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 20 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. .
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 where i = 1, . . . , M n − 1, δz i is the cumulative error in depth. At z 1 we have the error δz 1 = δ 1 . If δ i , i = 1, . . . , M n , are normally distributed random variables independent of each other and with the mean values and variances σ 2 z , δ i ∼ N ( , σ 2 z ), the relation between local (δ ) and cumulative (δz) depth errors is δz = Rδ .
Temperature perturbation value δτ n i can be expressed using the temperature gradient and the position error δz i as 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 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 E [δk] values are rather negative than positive.
In line with Eq. (18), the covariance in δk is 15 The variances of δk are found on the diagonal of Cov[δk] and tell how the variances σ 2 z are magnified through δ and to the variances in thermal conductivity estimates. 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 20 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 case L is assumed to be equal to 1 m, is 1.3 mm.

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 5 matrix B is defined as: or in matrix form: with the elements B n ij in B calculated and ordered in the same manner as in A in Eq. (10) and (12). A perturbation δρ can 10 be interpreted as a perturbation in the temperature δτ = −Bδρ in the first term of (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.

15
The behavior of the objective function F τ, (k, ρ) in (6) and (7) close to the optimum k * , ρ * is determined by the Hessian matrix H: where A, B, W and γ are defined in (11), (27), (7) and (6) and I is the identity matrix. A small perturbation δχ = (δk, δρ) of k * , ρ * will change F τ, by δχ T Hδχ. The Hessian H is positive definite with positive eigenvalues and orthogonal eigenvectors 20 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

25
The subsurface glacier profile observed in the cores consists of snow and firn with multiple ice lenses (Figure 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 ( Figure 1).
The evolution of subsurface temperature τ n i measured during five periods used for the simulations is shown in Figure 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 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 Figure 3B. The values consistently increase with depth at the rate of ca 0.11 W (m 2 K) −1 on average and somewhat slower in the fall of 2014. The temporal change in k * values can also be assessed with reference to Figure 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 year) −1 . Provided that the surface accumulation rates at the field site are 25 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 of 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 one standard deviation 30 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 ( Figure 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 increase of 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 5 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 Figure 4, the optimal k * (A) and ρ * (B) values are computed for the spring 2012 domain for different γ. The result varies for the extreme values of γ but is consistent for intermediate γ. The k * and ρ * values are indistinguishable for γ = 10, 10 3 and γ = 10, 10 3 , 10 7 , respectively. The eigenvector for the smallest eigenvalue ( Figure 4C) shows that the uncertainty is largest for 10 the deepest k * values when γ ≥ 10 in agreement with Figure 3B. When γ = 10 −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.

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 15 measurements τ n i are presented in Figure 5. The black markers close to the right border of each dataset show the positions of k * points. The color at any specific point in depth z i and time t n 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 n . 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 20 to constrain the optimization routine. For the first four calculation domains (spring 2012-2014 and fall 2014) this 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 color coding in Figure 5 by a factor of √ q = 3.

25
Optimized k * values are not very sensitive to single errors in subsurface temperature data: for j = 1 . . . 7 (panels A -G in Figure 5) the expected response of the thermal conductivity values k * j to a 1 • C error varies between −1.1·10 −3 and 1.1·10 −3 W (m K) −1 , corresponding to theĀ τ j − 3σ Aτj andĀ τ j + 3σ Aτj , whereĀ τ j and σ Aτj are the mean and standard deviation of the A τ j values. However, as it was demonstrated earlier (see the error bars in Figure 3B), a systematic time-independent bias in temperature data can result in notable deviations of the k * estimates. Several patterns in the spatiotemporal distribution of the sensitivities A τ can be noted. In most cases the sign of the sensitivity is positive between the depths of k * j−1 and k * j and negative between the depths of k * j and k * j+1 if k * j is the k value being tested for sensitivity (depth is marked by the black circle in Figure 5). This pattern is reversed when considering depth 5 levels further away from the circles (between k * j−2 and k * j−1 and between k * j+1 and k * j+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 of the temperature gradient just above and decrease just below that depth. These changes in temperature gradient are respectively compensated by negative and positive 10 deviations in k values. Due to the piecewise-linear interpolation of k profile based on J k * values, a perturbation in k * j also affects the thermal conductivities below and above it. Therefore k * j−1 and k * j+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 Figure 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 15 by the more intense colors in the sensitivity fields found in 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 colors around the circular markers in panels E -H of Figure 5 compared to panels A -D in the same figure.
4.4 Sensitivity of k * to errors in depths of temperature measurements 20 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 Figure 6, where different panels show results for the five simulation domains. The color 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 fall 2014) are significantly larger than for other k * j nodes due to 25 propagation of the corresponding high A τ values ( Figure 5 G, H) to the A z matrices as is shown in Eq. (15) and (21). These results are not shown.
In Figure 6 blue anomalies in 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 by lower k * values to keep the same heat flux. Secondly, the 30 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 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+1 and k * j−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 Figure 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 5 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 the larger A z values there.
Assuming that thermistor strings are in contact with the borehole walls every L = 1 m (see Eq. (25) Table 1 will increase almost by a factor of two. The results for the variance in Table 1 using Eq. (24) show how the variance σ 2 z in the positional error is magnified as a variance of the error in k. If σ z ≈ then Var[δk] ≈ 0.01 15 for most k-values, a rather small standard deviation.

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 Figure 7. All panels in the figure are dominated by colors 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 20 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 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 25 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 Figure 7 can be scaled by any assumed value E[δρ] to derive errors E[δk] to be expected from k * values. It follows that a bias of 50 kg m −3 will result in k * deviation of up to 0.1 W (m K) −1 .

Comparison of k * with earlier published results
The relation between optimized effective thermal conductivities and densities is shown in Figure 8 using markers and the associated linear fit: is illustrated by the thick black line. Also shown are the predictions from published similar fits. The Sturm et al. (1997) 5 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.

10
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)  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).

25
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 2 year −1 in the top 30 m of firn pack Greenland Summit. 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, page 400), the thermal diffusivity (κ * ) is calculated as κ * = k * (ρ * C) −1 . The resulting κ * values lie in the range from 7.55 to 30 48.77 m 2 year −1 with the mean value and standard deviations of 25.48 and 7.52 m 2 year −1 respectively, providing a prominent match with the results from Giese and Hawley (2015). It can be noted that quantification of k * using optimization technique requires a much less extensive dataset in terms of time and depth coverage. Furthermore, since infiltration of melt water 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.

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 con-5 ditions. 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, 10 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 result 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 15 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 climate similar to the one 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 be underestimated when simulated 20 using the parameterizations by Sturm et al. (1997) and Calonne et al. (2011).

Data and code availability
The data and computer codes used in the present manuscript are available from the authors upon request.  (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.

E[δk]
k * 1 k * 2 k * 3 k * 4 k * 5 k * 6 k * 7 k * 8     The color 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 Figure 2.    (Sturm et al., 1997), (Calonne et al., 2011) and .