1. Introduction

Nearly the entire Antarctic ice sheet is covered by a thick layer of firn (50–150 m), the transitional product between fresh snow and glacier ice (Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others, 2011). The ice underneath the firn layer has recently attracted most attention, due to its mass loss and subsequent contribution to sea-level rise (Reference StockerStocker and others, 2013). A common way to measure ice-sheet mass balance is to convert remotely sensed surface elevation changes (volume changes) into mass changes (e.g. Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009; Reference Helm, Humbert and MillerHelm and others, 2014). To convert from volume to mass, knowledge of spatio-temporal variations in firn layer depth/density is required (Reference Arthern and WinghamArthern and Wingham, 1998; Reference Ligtenberg, Kuipers Munneke and Van den BroekeLigtenberg and others, 2014). As in situ measurements in Antarctica are scarce, firn layer behavior is often simulated using a firn model (e.g. Reference HelsenHelsen and others, 2008). These models simulate firn layer processes (e.g. densification, vertical heat and liquid water transport) and are forced by the local climate at the surface. Because of the lack of firn densification observations, firn models are often poorly evaluated.

A number of firn cores across the Antarctic ice sheet provide valuable data to constrain the average depth and density of the firn layer (Reference Kaspers, Van der Wal, Van den Broeke, Schwander, Van Lipzig and BrenninkmeijerKaspers and others, 2004; Reference Van den BroekeVan den Broeke, 2008; Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others, 2011). These are point measurements in time, however, and do not contain information on rates of firn compaction. Surface elevation time series could potentially fill this gap, but they include several other signals, such as bedrock movement (Reference GunterGunter and others, 2014), ice-dynamical effects (Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009), surface melt, and accumulation variability (Reference Horwath, Legrésy, Rémy, Blarel and LemoineHorwath and others, 2012). Moreover, satellite altimeters either lack adequate temporal resolution (laser altimetry (e.g. Reference RivaRiva and others, 2009; Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012)) or contain significant uncertainties due to unknown surface penetration effects (radar altimetry (e.g. Reference Legrésy and RémyLegrésy and Rémy, 1998; Reference Ligtenberg, Horwath, Van den Broeke and LegrésyLigtenberg and others, 2012)). Borehole optical stratigraphy (Reference Hawley and WaddingtonHawley and Waddington, 2011), neutron probe observations (Reference Morris and WinghamMorris and Wingham, 2014) or fixing measurement devices in the firn column at various depths (Reference Hamilton, Whillans and MorganHamilton and others, 1998; Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others, 2010) provide direct measurements of firn compaction rates, but with limited spatial coverage. To assess the spatial and temporal scale of firn compaction rates, a method is needed that can be easily repeated without being destructive and can cover a large area.

In the companion paper (part I: Reference Medley, Ligtenberg, Joughin, Van den Broeke, Gogineni and NowickiMedley and others, 2015), such a method is described for obtaining firn compaction rates using airborne radar data from NASA’s Operation IceBridge mission over the Thwaites Glacier region, West Antarctica. Ground-penetrating radar data are often used to track isochrones (time horizons) in firn, which are then converted into accumulation rates (e.g. Reference VerfaillieVerfaillie and others, 2012; Reference MedleyMedley and others, 2013; Reference MiègeMiège and others, 2013). When such a radar is repeatedly flown over the same location, the temporal evolution of layer thicknesses can be converted into compaction rates (see Reference Medley, Ligtenberg, Joughin, Van den Broeke, Gogineni and NowickiMedley and others, 2015). In this paper (part II), these measured firn compaction rates are compared with model simulations. By doing so, we are able to examine firn compaction variability in three dimensions (temporal, vertical and along-track) and better constrain the conversion of volume changes into mass changes.

In Section 2, the firn densification model and its specifications are introduced, along with a brief description of the radar measurements. Thereafter, observed and modeled firn compaction rates are analyzed, compared and discussed. The paper concludes with a short summary of the main results.

2. Methods

2.1. Airborne radar measurements

We use measurements of firn compaction rates, observed by an airborne ultra-wideband microwave radar developed by the Center for Remote Sensing of Ice Sheets (CReSIS), University of Kansas, USA. The processing and specifications of the radar data are described in Reference Medley, Ligtenberg, Joughin, Van den Broeke, Gogineni and NowickiMedley and others (2015) and are only briefly discussed here. In a previous study (Reference MedleyMedley and others, 2013), annual horizons (isochrones) were tracked in the upper ~50 m of the firn layer of the Thwaites Glacier catchment. At locations with repeat-track measurements (fig. 1 of Reference Medley, Ligtenberg, Joughin, Van den Broeke, Gogineni and NowickiMedley and others, 2015), the evolution of the thickness between two tracked isochrones is converted into a firn compaction rate. Firn compaction is given as a positive value, i.e. a decrease in layer thickness. Because the uncertainties in isochrone depth would dominate over the relatively small compaction between two subsequent annual isochrones, thicknesses of 5 year layers (e.g. between isochrones dated to 1985 and 1990) were used to calculate compaction. In total, firn compaction rates are obtained for six 5 year layers, covering the period 1980– 2009/10. To reduce the noise signal in the radar data, while preserving its high horizontal resolution, the data were averaged along 1 km segments of the flight track. Radar-derived compaction rates are only used when >75 % (i.e. >750 m of the 1 km segment) were measured and found reliable (see Reference Medley, Ligtenberg, Joughin, Van den Broeke, Gogineni and NowickiMedley and others, 2015).

Fig 1. Observed (red) and modeled (blue) firn compaction rates, averaged over all 260 observations (a) and averaged over the time periods 2009–10 and 2010–11 (b). Grey and light-blue shading indicate the standard deviation in the radar observations (see Reference Medley, Ligtenberg, Joughin, Van den Broeke, Gogineni and NowickiMedley and others, 2015) and modeled compaction rates, respectively. Note that the observed and modeled average compaction rates in the upper ~5 m in (a) overlap as they are virtually identical.

2.2. Firn model

We use a one-dimensional, time-dependent firn densification model (FDM) to simulate the spatio-temporal evolution of the Antarctic firn layer. The model and simulation details are described in Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others (2011) and only briefly summarized here. The FDM calculates density, temperature, mass and liquid water content inside the firn layer, as well as their integrated effect on surface elevation. The thickness of a model layer in the FDM varies depending on the firn density, with a maximum of 0.2 m, while the total number of vertical model layers depends on the depth of the firn–ice transition (ρ = 917 kg m⁻³). At the surface, the model is forced with 6 hourly climate data from the regional atmospheric climate model RACMO2.1 (Reference Lenaerts, Van den Broeke, Van de Berg, Van Meijgaard and Kuipers MunnekeLenaerts and others, 2012a), including surface temperature, surface mass-balance (SMB) components (snowfall, sublimation, drifting-snow erosion/deposition and melt) and wind speed. For the region considered, the RACMO2.1 climate forcing is available at 27 km horizontal resolution and for the period 1979–2012. To obtain a realistic firn layer at the start of the simulation, a spin-up procedure is used that contains sufficient iterations of the 1979–2012 climate to refresh the complete firn layer (Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others, 2011). Firn compaction rates are stored for every model layer with a temporal resolution of 10 days.

The FDM is based on the principles of Reference Herron and LangwayHerron and Langway (1980), who divide firn densification into three stages. Near the surface (ρ < 550 kg m⁻³), the firn densification rate is largest, mainly due to packing of the (fresh) snow grains. Thereafter, ‘sintering’ is the main compaction process and the firn densification rate is roughly halved (Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others, 2010). When the density reaches the pore close-off depth (ρ = 830 kg m⁻³), densification only happens by compression of the existing air bubbles. To represent these three stages, we use the rate equations of Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010), tuned to fit 48 Antarctic firn density profiles, following Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others (2011):

(1)

where

(2)

(3)

where ḃ is the average annual accumulation (mm w.e. a⁻¹), is the average surface temperature over the entire period (1979–2012), g is the gravitational acceleration, ρ _i is the ice density (917 kg m⁻³), R is the gas constant, E_c and E_g are the activation energy constants (Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others, 2010) and T(z) is the local firn temperature at depth z. The density of fresh snow (ρ _s) is calculated following Reference HelsenHelsen and others (2008):

(4)

and depends on the local average annual surface temperature , accumulation ḃ (mm w.e. a⁻¹) and wind speed . The FDM also contains a firn hydrology scheme for the simulation of liquid water processes. However, surface melt does not occur in the region of the repeat radar tracks and is therefore not considered in this study.

In order to compare modeled and observed compaction rates, some assimilation of the FDM results is required. The FDM gridpoint closest to the observation in both space and time (18 October 2009, 19 November 2010 and 9 or 12 November 2011) is chosen, so no interpolation method is used. To obtain the total compaction rate between two isochrones (e.g. 2005 and 2000), the modeled densification rate is integrated between the radar-measured depths of these two layers. To mimic the fact that these layers move down with time, their depth is linearly interpolated between the start (in 2009 or 2010) and the end (in 2010 or 2011) of the observation period. Finally, all modeled 10 day compaction rates are summed over the period considered (2009–10, 2009–11 or 2010–11). This calculation is done for all six layers, which each encompass 5 years of accumulation (e.g. 1985–90), apart from the top layer, which represents the surface to the 2005 isochrone.

3. Results

As described in Section 2.1, the observed compaction rates are averaged in 1 km batches, resulting in 260 data points divided over three periods: 63 for 2009–10, 151 for 2009–11 and 46 for 2010–11. We find that the average compaction-rate profile derived from the 25 corresponding FDM gridpoints agrees very well with the average of all observations (Fig. 1a). The largest firn compaction rates are found near the surface for two reasons. First, the firn density near the surface is usually <550 kg m⁻³, resulting in higher compaction rates due to grain settling. Second, due to the nonlinear temperature dependence in Eqn (1), rapid compaction of fresh snow occurs near the surface in response to relatively high summer temperatures, compared to winter. In the second layer (6–10 m), the same two processes controlling compaction are still at work, but less pronounced. The seasonal variations in firn temperature decrease with increasing depth, and the transition to the second, slower stage of firn densification (at ρ = 550 kg m⁻³) occurs around 7–10 m in this part of West Antarctica (Reference Ligtenberg, Helsen and Van den BroekeLigtenberg and others, 2011). Both processes lead to a compaction rate that is roughly half that of the upper layer. Deeper down, firn compaction rates are smaller and more constant with depth. Below 10 m, hardly any seasonal temperature signal remains and variations in compaction rate are only caused by differences in local density (Eqn (1)). As density increases with depth, compaction rates decay to zero as the density approaches that of ice.

Figure 1b compares the modeled and observed compaction rates for two separate time periods (2009–10 and 2010–11). The smaller (larger) near-surface compaction rates over the 2009–10 (2010–11) interval are partly captured by the FDM. Within the survey area, both periods experienced similar temperatures (Fig. 2), while 2009–10 was wetter (428 mm w.e. a⁻¹) than 2010–11 (392 mm w.e. a⁻¹), thereby contradicting the finding that higher temperatures and accumulation rates lead to higher compaction rates. However, it is likely that the difference in the observed compaction rate in a certain period is also influenced by the climate conditions experienced immediately preceding that period, rather than climate conditions during the period itself. The year 2009 was indeed colder and drier than the long-term average, resulting in a relatively cold and dense firn layer at the start of the 2009 survey. As the radar surveys are conducted in late austral spring, the major part of the compaction takes place in the first part (summer) of the surveyed periods, which in the 2009–10 case results in lower compaction rates. Preceding the 2010–11 survey, the year 2010 was slightly warmer and received an average amount of winter accumulation, resulting in a warmer firn layer that contained more fresh, low-density snow at the start of the 2010 survey, compared to the firn layer at the start of the 2009 survey.

Fig. 2. Monthly surface mass balance (SMB, blue) and 2 m temperature (T_2m, red) for the period 2006–12, averaged over the survey area (similar to red area in Fig. 4a) from RACMO2.1. Firn air content (black) as simulated by the FDM, averaged over the same area. The orange lines indicate the timing of the radar surveys, and colored numbers show the October–October averages in SMB and T_2m for the corresponding year.

If correct, the above hypothesis would in part explain why the FDM is able to capture only part of the variability in compaction rates between 2009–10 and 2010–11. Two processes related to the relative amount of snow accumulation lead to temporal differences in firn compaction rate. First, in high-accumulation regions, more mass per unit time is added on top of the existing firn layer, thereby increasing the overburden pressure on the firn and enhancing the firn densification rate. Second, high accumulation results in quick burial of fresh snow, resulting in a thicker firn column that contains more firn to compact than thinner firn layers in low-accumulation regions. The firn densification rate in the FDM depends on the average annual accumulation (ḃ in Eqn (1)) rather than the overburden pressure. Therefore, the FDM is insensitive to temporal changes in loading (process 1). As a result, the variability in the modeled densification rate is represented only by the amount of low-density accumulation added (process 2), so total variability is underestimated. Moreover, snow microproperties (e.g. grain size) also play a role in the rate of compaction (Reference Freitag, Wilhelms and KipfstuhlFreitag and others, 2004), but are not included as prognostic variables in the current FDM. These properties, however, are likely to vary mainly between different snowfall events and are relatively constant on the seasonal to annual timescales considered here.

The standard deviation in the measured firn compaction rates (grey shading in Fig. 1a) is large, compared to the modeled one (blue shading). This finding is confirmed in Figure 3, where observed and modeled compaction rates are compared along a ~80 km transect (B in fig. 1 of Reference Medley, Ligtenberg, Joughin, Van den Broeke, Gogineni and NowickiMedley and others, 2015). The radar observations vary substantially on the kilometer scale, which could be due to either uncertainty in the radar measurements (0.05–0.1 m a⁻¹ for the upper two layers and 0.1–0.2 m a⁻¹ for the lowest layers) or a real signal forced by local climate variability not captured by RACMO2.1. Negative compaction rates (e.g. expansion) are observed at lower depths and deemed physically unrealistic. These are likely caused by the uncertainty in layer picking from the radar echograms. However, when averaged over the entire transect, the two methods agree well (table in Fig. 3). In the near-surface layer (orange), a substantial trend in observed and modeled compaction rates along the transect is found, with higher rates in the east. The eastern part of the transect is also located further south (see Fig. 4a), where higher accumulation rates and/or higher temperatures explain the higher compaction rate (Fig. 5). This trend is partly captured by the FDM, with higher compaction rates near the eastern end of the transect (around –104.1° E). In the other layers, no significant along-transect trend is found.

Fig. 3. Observed (dots) and modeled (lines) firn compaction rates along an IceBridge flight line (transect B in fig. 1 of Reference Medley, Ligtenberg, Joughin, Van den Broeke, Gogineni and NowickiMedley and others, 2015). Colors indicate the compaction rates for different firn layers: surface–2005 (orange), 2005–00 (red), 2000–1995 (light blue), 1995–90 (blue), 1990–85 (green) and 1985–80 (black). Values on the right show the average observed (Radar) and modeled (FDM) compaction rates (CR) for the entire transect. For reference, 0.3° longitude is equivalent to a horizontal distance of ~25 km.

Fig. 4. Densification rate of the upper 25 m of the firn layer, as modeled by the FDM (a), and for the periods 2009–11 (b), 2009–10 (c) and 2010–11 (d). The black dots in (a) mark the flight-line points for which the compaction rate is observed. The colored dots in (b–d) indicate the observed compaction rate for the corresponding time period.

Figure 4 shows the spatial variability in modeled and observed compaction rates of all layers combined (i.e. the upper ~25 m of the firn layer). Unfortunately, all observations are recorded in an area with similar modeled compaction rates (0.35–0.45 m a⁻¹), making it difficult to fully assess the spatial variability. The observed compaction rates again show more variability than the modeled rates, indicating that the model resolution is too coarse to capture all variability. Despite this finding, there is qualitative agreement in the spatial pattern: In both the observations and the FDM, the eastern track shows rates in the range 0.3–0.4 m a⁻¹, while the western track shows higher values (0.4–0.5 m a⁻¹) for both periods 2009–10 and 2009–11 (Fig. 4b and c). For the last time period (2010–11), only one flight track is available, but the observed compaction rates are in the same range (0.3–0.55 m a⁻¹) as the FDM simulation (0.35–0.45 m a⁻¹).

Figure 5a and b show that spatial differences in compaction rate are forced predominantly by the average annual accumulation, defined as the total precipitation minus drifting-snow erosion and total sublimation (the sum of surface sublimation and drifting-snow sublimation). The observations from the northwestern flight tracks (black dots) are situated in a relatively wet and warm region, explaining the (on average) higher firn compaction rates in that region. As mentioned earlier, high firn compaction rates are associated with high accumulation because of two processes: overburden pressure and the amount of low-density firn. The average surface temperature only has a second-order impact on the spatial pattern of the compaction rate, although it is difficult to fully isolate the effect as higher temperatures coincide with higher accumulation rates in West Antarctica (Fig. 5b and c).

Fig. 5. (a) Average compaction rate of the entire firn layer over the period 1979–2012, as modeled by the FDM. (b, c) Average accumulation (b) and average surface temperature (c) over the period 1979–2012, as simulated by RACMO2.1 (Reference Lenaerts, Van den Broeke, Van de Berg, Van Meijgaard and Kuipers MunnekeLenaerts and others, 2012a). The black dots in (a) mark the flight-line points for which compaction rate was measured.

To summarize, the observed and modeled firn compaction rates in this sector of West Antarctica compare well. Averaged over all data points, the vertical profile of firn compaction rate is similar: ~0.2 m a⁻¹ in the first 5 m and an order of magnitude smaller at greater depths. Both the FDM and radar observations show higher compaction rates in 2010–11, compared to 2009–10, suggesting that temporal variations in compaction rate are at least partly captured by the model. Large-scale spatial variability is also well simulated, with higher compaction rates occurring in the western locations, a result of the gradients in regional climate. This comparison provides confidence that the gridded FDM results contain a realistic representation of the average firn compaction rate on the model gridcell scale. The observed high-frequency variability in the radar measurements is likely caused by combination of small-scale variability in local climate that is not captured by RACMO2.1 and uncertainties in the radar measurements.

4. Discussion and Perspectives

Although the radar observations and FDM results qualitatively agree, the spatio-temporal variability and differences therein show a large degree of uncertainty. For the radar observations (see Reference Medley, Ligtenberg, Joughin, Van den Broeke, Gogineni and NowickiMedley and others, 2015), the uncertainty in the upper layers is of the same order of magnitude as the signal (0.1 m a⁻¹), while at greater depths (10–30 m) it is an order of magnitude larger. These large uncertainties are mainly due to (1) the isochrone-picking procedure, (2) a depth-cumulative error originating from the chosen firn density profile, and (3) the radar resolution. The large variability in the radar observations (Fig. 3) is in part a result of these measurement uncertainties. It is promising, however, that the average of several radar observations agrees quite well with the FDM, suggesting no systematic bias.

At the same time, it is highly likely that part of the compaction-rate variability represents real spatial differences initiated by small-scale variability in local climate. Gradients in bed topography result in subtle differences in elevation and slope of the ice-sheet surface. These surface irregularities impact the SMB as near-surface wind speed, sublimation and drifting-snow erosion/deposition differ on the upwind and downwind sides (Reference Lenaerts, Van den Broeke, Scarchilli and AgostaLenaerts and others, 2012b). Small-scale variability (1–10 km) in accumulation in this region is apparent in Reference MedleyMedley and others (2013). Along a flight track, dune-like patterns in surface elevation are likely to initiate a similar pattern in accumulation and firn compaction rates. In Figure 3, the upper-layer compaction rates (orange dots) between –104.5° E and –104.1° E indeed exhibit an oscillation; the amplitude of this oscillation (0.03–0.04 m a⁻¹ or 20–25%) is in the same range as SMB differences between up- and downwind sides of small topographic features, due to enhanced sublimation and/or drifting-snow erosion (Reference Frezzotti, Gandolfi and UrbiniFrezzotti and others, 2002). The climate forcing of RACMO2.1 has a gridcell size of 27 km, while these topographic features are an order of magnitude smaller and therefore not resolved in the model topography.

Modeled firn densification rates depend on local firn density and temperature, the average surface temperature, and the average annual accumulation (Eqn (1)). In reality, the rate at which firn compacts is also determined by the overburden pressure, meaning that periods with high accumulation would enhance densification, while during drier periods it would decrease with time. The FDM is not sensitive to these changes in loading and therefore probably underestimates temporal variability. However, the temporal resolution of the current observations (1 year) is not sufficient to test this hypothesis.

The density of fresh snow in the FDM is assumed constant in time and is parameterized as a function of the annual average accumulation, surface temperature and wind speed (Eqn (4)). This introduces an additional uncertainty in compaction rate, especially near the surface. In reality, the density of fresh snow is not constant (i.e. every snowflake is unique), meaning that with every accumulation event an error in the upper-layer density is introduced. As the firn compaction rate is relatively high for low-density snow, this error will quickly diminish with depth.

The FDM results also depend on the accuracy of the RACMO2.1 climate forcing. The simulated climate from RACMO2.1 is used to add (remove) mass to (from) the surface, but is also used for the calculations of the heat diffusion, fresh snow density and firn densification rate. RACMO2.1 has proven to realistically simulate the (near-) surface climate of the Antarctic ice sheet (Reference Lenaerts, Van den Broeke, Van de Berg, Van Meijgaard and Kuipers MunnekeLenaerts and others, 2012a; Reference Shepherd, Ivins and GeruoShepherd and others, 2012; Reference Sanz Rodrigo, Buchlin, Van Beeck, Lenaerts and Van den BroekeSanz Rodrigo and others, 2013), as well as that of separate regions (Reference Lenaerts, Van den Broeke, Scarchilli and AgostaLenaerts and others, 2012b; Reference DasDas and others, 2013; Reference Trusel, Frey, Das, Kuipers Munneke and Van den BroekeTrusel and others, 2013). In the Thwaites Glacier region, the annual average accumulation from RACMO2.1 agrees well with observations, but the interannual variability is less accurate (Reference MedleyMedley and others, 2013). As a consequence, while the total depth of the firn layer may be quite accurately simulated, the depth of individual annual layers might not be. As the observed compaction rates for different vertical layers are based on the dating of these annual layers, an error is introduced when the modeled firn compaction rates are integrated over these observed depths. This error appears small, based on the agreement in magnitude of the firn compaction rate when averaged for separate periods (Fig. 1b). However, if the temporal or vertical resolution of the radar surveys increases, this could become more serious.

Increasing the number of observations would greatly assist in firn model evaluation. Operation IceBridge only started flying in 2009, but improvement to the temporal resolution is not expected as the surveys are flown annually. Repeat flight tracks available from different climate regimes (i.e. SMB and temperature) in Antarctica and Greenland would certainly provide more insight into spatial compaction rate variations. Calculating compaction rates from repeat-track airborne radar, however, is time-consuming, as manual layer picking has to be checked thoroughly. There are, however, signs that (semi-)automatic layer pickers could substantially decrease the amount of work (e.g. Reference PantonPanton, 2014).

All these uncertainties aside, the results in this paper show that measuring firn compaction rates from repeat-track airborne radar allows an independent evaluation of firn model simulation and is, therefore, a very welcome addition to the field of firn modeling. Until now, it was not possible to verify modeled compaction rates over larger areas and for different time periods. In the part of West Antarctica discussed here, the observed and modeled compaction rates agree well when averaged to the model scale.