Firn cores

We select two test sites that are representative for the range of firn conditions found on the GrIS. At the first site, the North Greenland Eemian Ice Drilling (NEEM) site (Fig. 1a), a shallow firn/ice core was drilled in August 2007 to a depth of 80 m. The NEEM site is located in northwest Greenland, at 2484 m a.s.l., and experiences seasonal melt only in extremely warm years, such as the summer melt in 2012 (Reference NghiemNghiem and others, 2012). Based on the ice-core record, the annual mean accumulation rate is 22.7 cm ice equivalent and the mean temperature at the site is −28°C (Reference Steen-LarsenSteen-Larsen and others, 2011). The second site, the Flade Isblink ice cap, is located in northeastern Greenland (Fig. 1a), where an ice core was drilled in the summer of 2006. Located only 600 m a.s.l., this site experiences heavy seasonal melt. From this ice core, the mean annual accumulation has been estimated to be ∼50 cm ice equivalent (Reference LemarkLemark, 2010).

Fig. 1. (a) The elevation contours of the Greenland ice sheet and the location of the NEEM and Flade Isblink firn-core sites (marked with stars), along with the location of the EGIG line. (b, c) The radar data obtained with the ASIRAS instrument along the EGIG line in 2006 and 2008. The radar data cover ∼750 km in horizontal distance to a depth of ∼15 m. The data show the radar-return characteristics of the percolation zone (noisy area) and the dry-snow zone (low noise, stratified area) of the GrIS. Two distinct features are marked in both radar images; the features appear deeper by two layers in the 2008 image compared with the 2006 image. The image-derived firn-compaction rate from these two images is shown by the comparison of ten layers in each image. The signal-to-noise ratio is high on the image-derived compaction rate.

For both cores, density has been measured at 55 cm resolution and dated using high-resolution δ ¹⁸O measurements. The low resolution of the density measurements is not ideal for detailed studies of firn compaction, but since the δ ¹⁸O measurements provide a chronology with sub-annual resolution, it is sufficient to validate our methodology to assess firn-compaction rates inferred from ASIRAS data.

ASIRAS radar measurements

Figure 1b and c show the airborne Ku-band SAR altimeter ASIRAS data along the EGIG line (Fig. 1a). Before the ASIRAS data can be used to identify internal layers, radar echoes with an across-track mispointing above ±1° must be removed, as power peaks from the internal layers will be smeared out in these echoes. The remaining radar echoes are retracked individually to derive a surface, before a local-maximum algorithm is applied on each echo peak to obtain the depth relative to the derived surface. Next, the flight profile is divided into segments ∼1 km along-track, resulting in an average of 400 radar echoes per segment. For each segment, the number of detected peaks for 2 cm depth intervals is counted and normalized with the number of echoes in the segment. This gives the probability for detection of a peak in a given depth interval. Finally, Gaussian distributions are fitted to each layer-probability peak to obtain the average depth and standard deviation of all detected layers. From the ASIRAS data obtained along the EGIG line in 2006 and 2008, the standard deviations of the fitted Gaussian distributions are between 1 and 15 cm, with an average standard deviation of 7 cm in the dry-snow zone. In the percolation zone, the standard deviation of the fitted Gaussian distributions is between 3 cm and 1.2 m, with an average standard deviation of 15 cm (Reference StensengStenseng, 2011). In previous studies (Reference Hawley, Morris, Cullen, Nixdorf, Shepherd and WinghamHawley and others, 2006; Reference HelmHelm and others, 2007; Reference De la PeñaDe la Peña and others, 2010) a number of ASIRAS waveforms were averaged and peaks were identified manually. The data shown in Figure 1b and c are based on a new automatic method used to identify annual layers in the ASIRAS data from the firn (Reference StensengStenseng, 2011).

These layers are isochrones spaced 1 year apart, as seen when comparing the distinct features in the 2006 and 2008 images in Figure 1. The two features, marked with long arrows in Figure 1, appear two layers deeper in the 2008 image when compared with the 2006 image. Therefore, the layers inform us about the chronology of the uppermost firn in the dry-snow zone. The two images also provide a first-order estimate of the firn-compaction rate between measurements, as seen when the lengths of the two brackets, marking the same layers, are compared. This ASIRAS dataset provides new information about firn-compaction rates.

HIRHAM5 RCM output

In accordance with Reference SørensenSørensen and others (2011) the monthly mean output fields of the HIRHAM5 RCM (Reference Lucas-Picher, Wulff-Nielsen, Christensen, Aðalsgeirsdóttir, Mottram and SimonsenLucas-Picher and others, 2012) are used to force the firn-compaction model. The HIRHAM5 RCM is forced at the lateral boundaries with the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim re-analysis product (Reference Simmons, Uppala, Dee and KobayashiSimmons and others, 2007; Reference DeeDee and others, 2011) at T255 resolution (∼0.7° or ∼77 km) for the period 1989–2009. We performed a continuous RCM simulation with the HIRHAM5 RCM at 0.05° (∼5.55 km) resolution on a rotated grid, which is interpolated onto a 5 km equal-distance grid over Greenland. The HIRHAM5 RCM output is interpolated onto the geographical points of interest for the firn model using a linear nearest-neighbor interpolation scheme.

Firn compaction

The firn-compaction model used in this study is an extension of the firn model described briefly by Reference SørensenSørensen and others (2011), which aimed at better understanding measured elevation change from remote-sensing altimetry. A key assumption in the model is that meltwater is stored in the surface layer formed at the same time as melt occurred. Therefore, no percolation or latent heat release from refreezing is accounted for in that model. Here the firn-compaction model is updated in order to allow percolation of water into the snowpack. This allows for redistribution of mass and energy from the surface to internal layers in the firn column.

The classical description of firn compaction divides the firn into three stages (Reference Herron and LangwayHerron and Langway, 1980). In the first stage, defined down to the so-called critical density (ρ _c = 550 kg m⁻³; Reference BensonBenson, 1962), the firn-compaction process is assumed to be dominated by grain settling and packing. Below the critical density the compaction is believed to be controlled by power-law creep (Reference Arnaud, Barnola, Duval and HondohArnaud and others, 2000). The end of the second stage is not as clearly defined as the end of the first stage, where a clear transition in the firn-compaction rates is observed. In the literature, some authors define the third stage at the density range 804– 830 kg m⁻³, where the open pores in the firn become closed and the compaction is controlled by minimizing enclosed air volume. However, despite the change in physical properties, the third stage is not relevant to modeling the response to fast climate changes in order to estimate elevation change. Therefore, the compaction-rate factor of the second stage will be used from the critical density until the bottom of the firn column, at the ice density (ρ _i = 917 kg m⁻³). Based on the formulation given by Reference Herron and LangwayHerron and Langway (1980) the rate of compaction is modeled as

(1)

where c ₀ and c ₁ are stage-dependent compaction-rate constants. Reference Herron and LangwayHerron and Langway (1980) gave the compaction-rate constants with respect to water-equivalent accumulation rates. Here we follow the definition of compaction rates of Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010), where the accumulation rate, , is in units of kg m⁻² a⁻¹. In this framework, the compaction-rate constants based on Reference Herron and LangwayHerron and Langway (1980) are

(2)

where ρ _w is the density of water, R is the gas constant and T _a is the mean annual temperature at the site. The compaction-rate constants given by Reference Herron and LangwayHerron and Langway (1980) are based on observations of firn-core densities from the GrIS and the Antarctic ice sheet. The model has now been updated by several later studies to enable a time-dependent use of a firn-compaction model, that allows variable (temperature-dependent) firn-compaction rates (e.g. Reference Li and ZwallyLi and Zwally, 2004; Reference HelsenHelsen and others, 2008; Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others, 2010; Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others, 2011), but the original parameterization (Eqn (2)) still shows a good fit to sites on the GrIS, and the compaction rate is generally applicable for obtaining the average steady-state compaction rate. Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010) based their updated compaction-rate constants on a semi-empirical coupling of Nabarro–Herring creep (i.e. lattice diffusion; Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others, 2010) and normal grain growth, and tuned to fit three time series of direct measurements of firn-compaction rates from the Antarctic Peninsula. Here, the best-fitting expressions for the compaction-rate constants were found to be

(3)

where g is the gravity constant and T is the temperature at a given depth and time in the firn column (Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others, 2010). This more recent formulation of the compaction-rate constants can be used to derive a dynamic compaction model forced by changing the surface temperature and accumulation, but the parameters are based on observations from a limited area in Antarctica and may not be directly applicable to the GrIS. Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others (2011) used such a dynamic compaction model forced with temperature and accumulation rates from a RCM, and updated the compaction-rate parameters in Eqn (3) using a large set of Antarctic firn-density profiles. Their updated best-fit parameters are, however, influenced by the assumed climate forcing from the RCM, which may differ from the in situ climate, and therefore not be generally applicable. We use a similar approach to that of Reference SørensenSørensen and others (2011) and Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others (2011), and develop a dynamic compaction model forced by a RCM. We use the model parameters of Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010) as our starting point, as they are based directly on in situ climate observations, and, as explained below, we take into account possible deviations of the RCM climate forcing field from the in situ conditions in the tuning of model parameters to density profiles. Here the accumulation rate, , is determined at each gridpoint in the model by the weight of the overlying snow/firn and the time since burial. As mentioned above, the Reference Herron and LangwayHerron and Langway (1980) compaction-rate constants can simulate the observed firn-density profiles from dry-firn cores on the GrIS and, as the seasonal oscillation of the surface temperature gets damped with depth, the firn compaction of the second stage should approach isothermal conditions, as modeled by Reference Herron and LangwayHerron and Langway (1980). The difference between the two formulations as T → T _a is given by

(4)

If this expression for γ is added to the second stage of Eqn (3), the Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010) formulation will (in line with the Reference Herron and LangwayHerron and Langway (1980) formulation) have an accumulation-dependent description of firn compaction under steady-state conditions. Equation (4) also indicates that the Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010) formulation of firn-compaction rate compacts the firn too fast compared with observed conditions in Greenland. Since the Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010) compaction-rate factors are based on measurements of dynamic firn compaction, rather than being derived from static observation from ice cores (e.g. the Reference Herron and LangwayHerron and Langway (1980) formulation), we use Eqn (3) in this paper, despite the difference between the Reference Herron and LangwayHerron and Langway (1980) and Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010) compaction rates for the second stage.

A number of empirical relationships for the surface density can be found in the literature for both Antarctic (Reference Kaspers, Van de Wal, Van den Broeke, Schwander, Van Lipzig and BrenninkmeijerKaspers and others, 2004; Reference HelsenHelsen and others, 2008; Reference LenaertsLenaerts and others, 2012) and Greenland (Reference Reeh and FisherReeh and others, 2005) surface snow. The Greenland study gives empirical parameterizations based on a compilation of study sites with 10 m temperature, T _f, ranging from −31 to −5°C,

(5)

where is the mean annual density. To limit the number of parameters in the firn-compaction model, we adopt this formulation by replacing T _f with the RCM modeled monthly mean surface temperature at the site.

Meltwater retention

Surface melt is significant on the GrIS (Reference Hall, Nghiem, Schaaf, DiGirolamo and NeumannHall and others, 2009), and when modeling firn compaction it is essential to include a description of refreezing of meltwater, in order to capture the dynamic range of firn compaction of the GrIS. The ice-sheet surface-energy balance is described by the RCM, therefore the retention scheme has only to account for the transport of liquid water present at the surface. A number of retention parameterizations can be found in the literature (e.g. Reference Coléou and LesaffreColéou and Lesaffre, 1998; Reference Janssens and HuybrechtsJanssens and Huybrechts, 2000; Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others, 2011). Here we account for the amount of liquid water distributed to the firn layers according to the maximum retention potential of each layer, R _max:

(6)

where R _energy is the maximum refreezing potential of a firn layer with initial temperature T and mass m, before it reaches its melting point, and R _pore is the maximum pore space in the firn. Given the specific heat capacity of ice, c, and the latent heat of fusion, L, the maximum refreezing potential can be defined as

(7)

Prior to formation of ice lenses, meltwater is accumulated and refrozen in the open pores of the firn. R _pore is defined by simple volume considerations:

(8)

where λ is the thickness of the layer and ρ _ip is the density of the firn where it becomes impenetrable for percolating water. Since limited information about the impenetrable density can be found in the literature, we perform a series of experiments (I1–I3; Table 1) to determine the value of ρ _ip most suitable for the firn-compaction formulation presented here.

Table 1. Inversion experimnents performed for the two firn-core sites and summary of the applied assumptions for each of the experiments

Assuming that meltwater will occupy the open pore space in the firn, the thickness of an individual firn layer is not affected by the percolation, but refreezing of meltwater will add mass to the layer and densify the layer. The density of a layer is then given by

(9)

where ρ_λ is the density of a layer with thickness λ, and RF is the amount of refrozen water. The exact location of ice lenses within the firn layer is not determined in this formulation.

As long as water is available and the porosity, η (where η = 1 – ρ/ρ _i), of the firn allows percolation, water that does not refreeze will percolate to deeper layers by gravitation. However, a small fraction of the water is able to resist gravity-induced percolation by capillary and adhesive forces, which is often referred to as irreducible water (Reference Coléou and LesaffreColéou and Lesaffre, 1998). The ratio of irreducible water, θ _irr, has been shown to have an exponential relationship with the porosity and is given by (Reference Schneider and JanssonSchneider and Jansson, 2004)

(10)

This results in an irreducible water mass given by

(11)

where is the irreducible water content (liquid water) left from the last percolation event and m _per is the percolating water entering the layer. Note that ηλρ _w is the maximum water content in a layer of thickness λ at a given time.

Temperature evolution

The derivation of the general heat equation is given by Reference Cuffey and PatersonCuffey and Paterson (2010), and assuming one-dimensional (1-D) heat transfer, it is written

(12)

where k _T is the thermal conductivity and f is internal heat production. We further assume (1) internal heat production (besides the latent heat release of refreezing, which is treated separately) is negligible and (2) the vertical heat transport is governed by the addition of layers at the top of the firn column. Then the heat equation (Eqn (12)) can be solved efficiently by numerical methods (Reference Schwander, Sowers, Barnola, Blunier, Fuchs and MalaizéSchwander and others, 1997). The above formulation of firn compaction is highly temperature-dependent (Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others, 2010), and temperature calculations using the general heat equation are density-dependent. The numerical implementation of the model solves the equations iteratively, until the solution converges.

Latent heat release from refreezing heats the firn layers. As the percolation of meltwater is assumed to be instantaneous, the heat release is also assumed to be instantaneous. Therefore, when water percolates through the firn column, the heat sources can be computed by updating the temperature profile between time-steps. The latent heat release of refreezing is estimated by the analytical solution for ΔT given by

(13)

This assumption makes the code less computationally demanding than solving the full heat equation, and it is accurate for the 1-D case that we use here.

Firn model initialization

Initialization of the dynamical firn model is needed, as the firn model describes the addition of firn layers on top of a pre-existing firn surface. As initial conditions, we use the steady-state firn-compaction model of Reference Herron and LangwayHerron and Langway (1980) with the inclusion of ice lenses (Reference ReehReeh, 2008) and the isothermal compaction-rate constants of Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010). The model assumes the 1989–2009 average climate from the HIRHAM5 RCM, and provides a mean density profile for the site for equally spaced layers in the firn. It is thereby a computationally efficient representation. With compaction rates given by the model, the steady-state density profile is interpolated onto the time domain used by the dynamical firn-compaction model.

The dynamical model starts in January 1989 (the start of the HIRHAM5 RCM record), by adding monthly layers to the surface. The combined initial steady-state firn model and dynamical firn model will affect the firn-temperature profile towards the annual mean for the first few months after January 1989, but the effect gradually deceases.

NEEM site: dry case

Experiment F1 is denoted the control experiment, since it applies all available information for the inversion, whereas F2–F6 evaluate the bias when limiting the model parameter space or the a priori information. Reducing the model parameter space by one parameter at a time (F2 and F3) results in a higher root-mean-square error (rmse) than in the control experiment, as shown in Table 3. This indicates that all the tuning parameters are necessary for the model to have the right degrees of freedom. Experiment F4 neglects the seasonality of the external forcing tuning parameters, which again results in an increased misfit to the observations. Reducing the number of model parameters to an absolute minimum (experiment F5) reduces the misfit compared with F1. However, the limited surface melt at the NEEM site results in a less constrained value of γ_{lq ,} which may increase the misfit. Therefore, when assuming γ_lq = 1, as in experiment F5, the unconstrained γ_lq does not affect the final result. Assuming γ_lq = 1 may be a good approximation at NEEM, but may not be valid for the entire GrIS firn-covered area. Finally, the F6 rmse is lower than found for the F1 experiment, which suggests that the initial inversion set-up (F1) is over-constrained, and that the information used in F6 is sufficient to constrain the inversion. The γ_{ρ s}, γ_{c 0} and γ_{c 1} model parameters for F6 are all well determined by the inversion and exhibit a low coefficient of variation (σ/γ ∼ 0.08), while the external forcing tuning parameters result in a higher coefficient of variation (σ/γ ∼ 0.22), which may be considered the accuracy of the inverted parameters. Further, the RCM output field provides a mean annual accumulation at the site of 18.4 cm, which is significantly lower than the annual accumulation of 22.7 cm ice equivalent from ice-core measurements. After tuning, the mean annual accumulation is 21.9, 22.1, 25.7 and 23.0cm for F1, F6, L1 and L2, respectively. For all the experiments except L1, this is within the error of the observation. In addition to the RCM output accumulation rate field being too dry, the results of the inversion indicate that the deviation of the RCM output accumulation rate is stronger in cold months than in warmer months at the NEEM site.

Table 3. The result of the inversion experiments performed at the three sites. The results for the ASIRAS radar inversion are only expressed as the mean of Figure 6b. For the EGIG-line ASIRAS result, the root-mean-square error (rmse) is derived as the mean of the 20 individual inversions performed in the test area, in units of air-equivalent radar travel (m)

The two additional inversion experiments, L1 and L2, are performed to investigate the bias of limiting the observations used in the inversion. Additional a priori information is needed to constrain the tuning of the second-stage firn compaction, γ_{c 1,} if the observational dataset is limited to densities less than ρ _c. The results of γ_{c 1} in F1–F6 all agree with the predicted value based on Eqn (4). Hence, a Gaussian a priori for γ_{c 1} is added to constrain the second-stage firn compaction based on Eqn (4) (Table 2).

Figure 3 shows the results for F6, L1 and L2. A common feature of the modeled densities is the deviation from the observations at 1 m below the January 1989 layer (marked by a horizontal line in the right panel). This is an artifact of the transition from the steady-state initialization to the dynamic firn-compaction model, which can be expected due to the model set-up. Both F6 and L1 result in similar densities in the top firn, while L2 tends to overestimate the surface density at the site and underestimate the compaction rate (Fig. 3; Table 3). This is also reflected in the increased coefficients of variance, σ/γ, of γ_ρs and γ_{c 0} for L2, which are almost double that of F6.

Fig. 3. The NEEM inversion results for experiments F6, L1 and L2, compared with the observations at the site (in red). The panels on the left show the symmetrical correlation matrices for each of the three experiments. The difference between the experiments is illustrated by the correlation matrices. The panel on the right shows both modeled densities and ages for the three experiments. The transition from dynamic to initialized firn modeling at the beginning of January 1989 is marked with the horizontal line, and ρ _c is marked with the vertical line.

The difference between the F6, L1 and L2 results for the NEEM site firn core can also be seen in the correlation matrices in Figure 3. The fully constrained F6 experiment shows a strong (< −0.8) negative correlation between the modeled parameter tuning for warm and cold months in accumulation, and for the surface density and c ₀, as illustrated in the upper left panel of Figure 3. This correlation indicates a trade-off between γ_{b w} and γ_{b c} and also γ_{p s} and γ_{c 0} in the accepted solution of the random walk. The correlations between the two accumulation parameters are constrained by the chronology, whereas the correlation between γ_{p s} and γ_{c 0} is constrained by the density information. This also explains why the L1 experiment only finds an anticorrelation between γ_{p s} and γ_{c 0} and L2 only finds an anticorrelation between γ_{b w} and γ_{b c,} as seen in the left column of Figure 3. As seen in Table 3, the L1 inversion obtains an optimal fit with observed densities by increasing the annual accumulation by ∼10% compared with F6, which results in a rmse of 1 year in the modeled age. The L2 experiment finds the optimal parameter set with an increased surface density, which introduces a higher rmse for density.

The rmse normalized with respect to the F6 rmse, gives L2 the lowest combined error followed by F6 and F5, while L1 has the highest normalized error, mainly attributed to the misfit of 1 year in the model age, as seen in Figure 3. The normalized rmse for L2 is smaller than for L1, which indicates a better constraint of the combined product of firn compaction and external forcing by only applying information about the chronology in the inversion, as in experiment L2.

Flade Isblink: wet case

The model results, shown in Figure 4, indicate the model is capable of reproducing the mean densities and the chronology at the Flade Isblink site in all experiments, but neither high- nor low-density peaks in the density profile are captured by the model. Both the minimum and maximum peaks in density can be explained as resulting from neglecting horizontal movement of percolating meltwater in the model and the time discretization of the model. Neglecting the horizontal movement of water flow affects the modeling result in two steps. First, horizontal water flow would be directed along ice lenses, increasing the density around these. Second, due to the missing addition of mass with depth, the inversion tries to find the optimal model parameters in order to fit the density profile by increasing the surface density (γ_{ρ s} = 1.47 for I3). Both steps result in a limited precision of the error estimate based on the chronology measure. In addition, the rather coarse time discretization of the model will limit the formation of high-density layers in the firn profile (as short-lived percolation events (e.g. rain or strong melt events) are not included in the model) which might form ice layers in the top firn. These layers may be additionally enhanced at later percolation events and thereby be an important feature in the formation of high-density layers.

The I3 experiment (ρ _ip = ρ _i) has the lowest density misfit for the inversions and is therefore favored as the best representation of the Flade Isblink firn column. This solution also shows two promising features: (1) A decrease of the modeled densities dated to the 1990s, which is also seen in the observations. This is a prominent feature that might be enhanced if horizontal water flow were allowed. (2) The correlation matrix of I3 shows an anticorrelation for the parameter pairs (γ_{b w}, γ_{b c}) and (γ_{lq w}, γ_{lq c}) and a less-noisy correlation matrix than I2 (Fig. 4). When the experiments are compared, ice lenses are the only real obstacle for the water percolating in the firn. Therefore, the initialization of the firn does not bias the dynamically driven layer above. As observed in Figure 4, the uppermost density is increased from experiment I1 to I3; this can be attributed to the run-off criteria in the model description, which removes the water when downward motion is ceased. The latent heat release at sites like Flade Isblink is an important heat source in the modeled temperature profile. A comparison with borehole temperatures at the site showed that the modeled temperature profile is biased towards lower temperatures, due to difficulties in estimating the cold content of the initialized temperature profile. This suggests that the temperature inversion is less constrained by the in situ parameters, which is expressed by the high standard deviation of the temperature tuning, compared with its a priori information

Fig. 4. The results of the Flade Isblink inversion experiments and observations at the site. The correlation matrices for each experiment are shown in the left panels. The right panel shows the modeled densities and ages.

The result for the RCM tuning indicates an overestimation of the precipitation by the RCM output files (Table 3). This suggests that the RCM tends to deliver too much precipitation at the margin of the GrIS (as seen at the Flade Isblink site), and be too dry in the interior, as observed with the NEEM site inversions. The high γ_{ρ s} combined with high temperatures at Flade Isblink results in high-density layers being deposited at the surface by the model. Since the mean density profile is captured by the model (Fig. 4) the high-density surface layers cannot explain the γ_{c 1} < 1 found in all of the experiments.

EGIG line: spatial case

Figure 5 shows the 2008 ASIRAS radar data along with the modeled depth of annual layers using the tuning parameters determined in the NEEM F6 experiment. The strong reflections in the radar data have previously been associated with deep hoar formation in the fall (Reference De la PeñaDe la Peña and others, 2010). However, such microphysical processes are not modeled here, and we assume that the annual stratification is represented by the modeled October layer. The modeled thickness of the annual layers are in good agreement with the first few observed layers, but the agreement tends to deviate with depth. In particular, the thinner layers on the eastern side of the GrIS ice divide (located at ∼35° W) are not captured by the modeled layers when using the parameters inferred from the NEEM F6 experiment. The deviations between the observed firn layers and the layers modeled using the parameters found for the NEEM F6 experiment emphasize the need for a separate investigation of the tuning parameters at the EGIG line. In the shaded box in Figure 5, the difference between the observed and modeled layers changes sign at ∼38° W. Several causes can be given for this (e.g. the HIRHAM5 RCM forcing could be biased, the firn-compaction model could be biased in the physical description or the spatial/temporal patterns of the tuning parameters could vary over the GrIS). To investigate this, an inversion for the tuning parameters is performed at each of the 20 HIRHAM5 RCM gridpoints in the selected area on the EGIG line. The inversion follows the method developed for the NEEM L2 experiment with ρ _ip = ρ _i, since the radar data only provide countable layers and no information on densities. In the selected area of the EGIG line, the slope of the isochrones is small and the seasonal melt is almost absent, reducing the uncertainty of the inversion. Figure 6a shows the result of the inversion for the most likely solution, along with the traced ASIRAS layers. All three Gaussian distributions found by the EM algorithm are shown as contours in Figure 6b, where the Gaussian distributions have been normalized in accordance with their probabilities found by the EM algorithm. The solution illustrated in Figure 6a is assigned a >50% probability by the EM algorithm for the majority of the inversion points. However, the two remaining solutions (not shown in Fig. 6a) at each point cannot be discarded as bad representations of the firn, with the biases found by the NEEM L2 experiment in mind. Quantitative comparison with Reference Morris and WinghamMorris and Wingham (2011, their fig. 2) shows that the three solutions of the EM algorithm capture the observed densities, with a slightly higher density in the high-density layers. This suggests that studies of surface density should be conducted at locations with a good instrumental record of forcing parameters, in order to be sure the inversion is well constrained and allowing the surface density to be the primary unknown in the inversion. However, this is beyond the scope of the present study.

Fig. 5. The 2008 ASIRAS radar image with modeled October layers estimated using the model parameters determined in the NEEM F6 experiment. The observations and modeled layers are in good agreement for the first few layers and then gradually deviate with depth. The thinner layers on the eastern side of the ice divide, especially, are not well represented by the model. The ice divide can be identified in the image by the change in slope of the layers at ∼35° W. The shaded area (∼37−39° W) was selected for further inversion studies (Fig. 6) using the ASIRAS radar data as observational input for the inversion method developed for NEEM site L2. The depth scale is in air-equivalent radar travel distance.

Fig. 6. The result of firn-compaction inversion using the ASIRAS radar data for the section of the EGIG line shaded in Figure 5. The mean rmse between the modeled and observed layers is only 9 cm, which is within the error of the layer tracing. (a) In blue are the traced layers from the 2008 radar image, and stars indicate the locations of the modeled annual layers at each of the 20 HIRHAM5 gridpoints. The depth scale is in air-equivalent radar travel distance. (b) Contour plot of the normalized probability distribution based on the EM algorithm for each of the 20 gridpoints. The grid numbers are shown in (a).

Table 3 lists the mean of the preferred tuning parameters of the individual inversions performed on the EGIG line. The rmse of the observed layer depth is 9 cm, which is equivalent to the standard deviation of the traced layers in the area. As seen in Figure 6a, the interannual variability in layer thickness along the 67.5 km of the EGIG line is captured by the model. This interannual variability in the layer thickness is not tuned by the inversion scheme, since the same seasonal tuning is used for all layers at a HIRHAM5 gridpoint. Hence, this suggests that the HIRHAM5 RCM output fields provide accurate temporal variation in the forcing field to the firn model. The average of the inverted parameters for the HIRHAM5 RCM at the EGIG line is similar to the NEEM results, both in absolute magnitude and in uncertainty, again indicating that the RCM output fields underestimate the inland precipitation on the GrIS.

The results of tuning the firn parameters, γ_{ρ s}, γ_{c 0} and γ_{c 2}, are not as straightforward to interpret. The limited uncertainty found for γ_{c 1} on the EGIG line may be attributed to the a priori information applied, and, as illustrated by the spread of EM solutions for γ_{c 0} in Figure 6b, the tuning parameters are not well determined and show high spatial variability. The average of the 20 individual inversion results gives γ_{c 0} = 1.4 and γ_{ρ s} = 0.73. The value of γ_{c 0} is therefore significantly higher than found at the NEEM site. The NEEM L2 experiment tended to overestimate γ_{v s} when only the chronology was applied as prior information (Table 3). This is not seen in the mean solution of the EGIG inversions, where the surface densities have to be increased to be comparable with the NEEM result. If the anticorrelation of the two parameters as found at NEEM is interpreted as a linear relationship, a 0.1 increase in γ_{ρ s} would only result in a reduction of γ_{c 0} to 1.27. Therefore, the lack of correlation between the two parameters cannot explain the significantly higher firn-compaction rates found on the EGIG line than at NEEM. This indicates a site dependency of the compaction-rate factors, which cannot be reproduced by the firn parameterization used here.