Modeling enhanced firn densification due to strain softening

. In the accumulation zone of glaciers and ice sheets snow is transformed into glacial ice by ﬁrn densiﬁcation. Classically, this process is assumed to solely depend on temperature and overburden pressure, which is controlled by the accumulation rate. However, exceptionally thin ﬁrn layers have been observed in the high-strain shear margins of ice streams. Previously, it has been proposed that this ﬁrn thinning can be explained by an enhancement of ﬁrn den-siﬁcation due to the effect of strain softening inherent to power-law creep. This hypothesis has not been validated, and the greater ﬁrn densities in the presence of horizontal strain rates have not yet been reproduced by models. Here, we develop a model that corrects the ﬁrn densiﬁcation rate predicted by classical, climate-forced models for the effect of strain softening. With the model it is conﬁrmed that strain softening dominates the ﬁrn densiﬁcation process when high strain rates are present. Firn densities along a cross section of the Northeast Greenland Ice Stream (NEGIS) are reproduced with good agreement, validating the accuracy of the developed model. Finally, it is shown that strain softening has signiﬁcant implications for ice core dating and that it considerably affects the


Introduction
Firn densification refers to the transformation of snow into glacial ice, which occurs in the uppermost layers of ice sheets 15 and glaciers within their accumulation zones, when old snow, respectively firn, is buried under younger snow. The overburden pressure gradually increases and causes a densification of the firn. Large scale ice flow is not considered by firn models even though it is known that ice is a non-newtonian material where strain reduces the viscosity (Goldsby and Kohlstedt, 2001). In this paper we demonstrate that that this effect can have a substantial impact on firn densification.
Firn densification is conventionally divided into stages where different physical mechanisms dominate. Initially, the Newto-20 nian grain-boundary sliding is dominant for densities of ρ ≤ 550 kg m −3 , which is referred to as stage 1 of firn densification (Alley, 1987). In stage 2, at densities of 550 kg m −3 ≤ ρ ≤ 830 kg m −3 , the non-Newtonian dislocation creep, also known as power-law creep, dominates the densification process until bubble close-off (BCO) (Maeno and Ebinuma, 1983). Beyond this point, defined as the firn-ice transition, the further densification is slowed as the enclosed gas-bubbles get compressed and eventually diffuse into the ice matrix (Salamatin et al., 1997). 25 softening. Based on a constitutive equation by Duva and Crow (1994) describing the behavior of porous media under powerlaw creep, Gagliardini and Meyssonnier (1997) have developed a flow model for alpine glaciers, which inherently considers compaction, and in particular strain softening, of the firn. This approach requires knowledge of the internal, horizontal stresses from flow modeling and is widely used in studies of alpine glaciers (Lüthi and Funk, 2000;Zwinger et al., 2007;Licciulli et al., 2020), but has not been applied to polar ice sheets. 70 In this paper we aim to model the effect of strain softening in the firn. In order to do so, a correction factor is derived that can be applied to correct the firn densification rate, predicted by any climate-forced firn densification model, for the impact of horizontal strain rates.

Theory
A firn densification model in a Lagrangian formulation expresses the densification rate of a firn layer as a function of external 75 forcing parameters and internal parameters, representing its current state. The external parameters are generally time-variable.
In a climate-forced model these are the temperature T and accumulation rateḃ, respectively the thereof derived overburden load σ. As an internal parameter most firn models only consider the current density ρ of the firn layer. Newly formed snow layers at the surface further require the initial snow density ρ 0 as a boundary condition. For a specific site this is however in most applications assumed to be constant with time. The densification rate of a climate-forced model (denoted with the subscript c) 80 is hence given in the form of with the time t. Following Morris et al. (2017) this densification rate relates to a corresponding volumetric strain rateε. As this type of model only considers densification in vertical direction the volumetric strain rate matches the vertical componentε zz,c , which hence reads as Because no additional horizontal strain rates are considered to affect firn densification in a climate-forced model, the strain rate tensor of such a model at a specific point within the firn is expressed bẏ strain rate tensor becomeṡ where the horizontal components (ε xx ,ε xy ,ε yy ) generally affect the vertical strain rateε zz by the effects of horizontal divergence and strain softening, such thatε zz =ε zz,c . More particularly this means that the corresponding densification rate is in consequence also altered by the horizontal strain rate components. Inserting Eq. 2 into Eq. 5 gives Accordingly, the densification rate with strain softening can be computed from the densification rate predicted by a climateforced model when a scale factor given by the ratio of the corresponding vertical strain rate components is known. We aim in 100 the following to determine this scale factor.

The scale factor
In firn stage 2, where power-law creep is dominant, firn densification can be modeled by the constitutive relationε zz = Aσ n with the creep factor A, representing the temperature impact, and the creep exponent n. The power-law creep equation can be generalized analogously to the derivation of Nye's generalization of Glen's flow law (Nye, 1957;Greve and Blatter, 2009, 105 Ch. 4), with the only difference being that stress is used instead of the stress deviator. This adaptation is valid without limitations and avoids the assumption of incompressibility that would be implied by the use of the stress deviator. The additional assumption of isotropy made in Nye's derivation also can be followed as firn is generally isotropic. The generalized power-law creep equation then reads aṡ 110 where the vertical strain rate of the densification is proportional to the causative overburden load, scaled by an effective viscosity η, given by This effective viscosity again depends on the effective strain rateε E = 1 2 ε 2 xx +ε 2 yy +ε 2 zz +ε 2 xy 1/2 with the exponent m = 1 − 1/n. These equations illustrate, how the horizontal strain rates lead to a decrease of the effective viscosity, which again for 115 the same load gives rise to a higher vertical strain rate. This is exactly what the term strain softening denotes. As the enhanced vertical strain rate is dependent on the effective viscosity which again depends on the vertical strain rate itself, its derivation is not straight forward. Eqs. 7 and 8 can be formulated similarly for the case of no horizontal strain. We note that the load and the creep factors are not affected by the horizontal strain rates and therefore must cancel out when we take the ratio in eq. 6 to obtain the scale factor.

120
The effective viscosity of a purely climate-forced model η c will be overestimated, as it neglects the effect of strain softening by horizontal strain rates. The effective strain rate according to the reduced strain rate tensor (3) isε E,c = 1 2ε 2 zz,c 1/2 , which leads to a reduced vertical strain rate oḟ By dividing Eq. 7 with Eq. 9, we obtain the scale factor and note that it is given by the ratio between the effective viscosity of 125 the climate forced model and the 'true' effective viscosity: Inserting the effective viscosities according to Eq. 8 with the corresponding effective strain rates giveṡ which can be rewritten in the form of whereby the following two variables are defined: While all components of r H are known, which are the horizontal strain rates derived from velocity fields and the vertical strain 135 rate of the climate-based model, the variable r V corresponds exactly to the scale factor that is sought. Hence, Eq. 12 needs to be solved for r V in order to obtain the scale factor for correcting the original densification rate of a climate-forced model for the effect of strain softening.
The solution of Eq. 12 depends on the creep exponent. Dislocation creep is the key process driving densification in the second firn stage (Maeno and Ebinuma, 1983). In this paper we therefore use a creep exponent of n = 4, characteristic for 140 dislocation creep (Goldsby and Kohlstedt, 2001;Bons et al., 2018). At this point only the solution shall be given, because its derivation would exceed the scope of the paper. We find, This scale factor is the strain softening enhancement. An alternative, simpler solution for n = 3 is derived in (Oraschewski, 2020).

Regularization
The densification in climate-forced models, such as the HL model, is characterized by an exponential decay towards the density of solid ice. The vertical strain rateε zz,c in such models goes to zero as the density approaches 917 kg m −3 . This means that r H , and thus the scale factor goes towards infinity. This singularity causes an unphysical behavior of the strain softening model, where the density of ice is approached almost instantaneously at a certain point.

150
To circumvent this issue, a regularization is introduced. Inspired by regularized Glen's flow law (Greve and Blatter, 2009, Ch. 4) a residual strain rateε 0 is added to the vertical strain rate to ensure a finite correction factor: Leaving the perspective of firn densification modeling, the residual strain rate can be associated with the general thinning of firn and ice layers in ice sheets that is induced by the flow of ice towards the ice sheet margins. While the densification part of 155 the vertical strain rates approaches zero, this contribution remains finite. The corresponding vertical strain rate is in the order ofε 0 ≈ −2 × 10 −4 yr −1 , which we choose as the residual strain rate in the following. The derivation of this residual strain rate can be found in (Oraschewski, 2020).

Implementation
The strain softening scale model is implemented as an optional module into the CFM by Stevens et al. (2020). Supported by its 160 Lagrangian modelling scheme, the implementation itself is rather simple and computationally cheap. The uncorrected vertical strain rate is computed from the initial densification rate according to Eq. 2 using a classical climate-forced densification model.
Together with the input data for the horizontal strain rates the variable r H can then be computed by Eq. 16, from which the correction factor of the vertical strain rate is computed with Eq. 15.
horizontal strain rates are loaded in the form of the principal horizontal strain rates (see e.g. Nye, 1959). In this way the shear components disappear, such that the amount of input data is reduced without the loss off any relevant information.
The strain-softening model is only applied in the second stage of firn densication for ρ > 550 kg m −3 as power-law creep is thought to dominate firn densification only in this range, while the Newtonian grain-boundary sliding is driving the densification process before. Although a smooth transition between the two processes over a range of densities is expected (Hörhold et al.,170 2011), its exact course is unknown and a sharp transition is assumed in our present implementation. Nonetheless, attempts to implement a smooth transition did not affect the model output significantly (not shown).
Within the CFM framework the strain softening model can be executed in combination with any of the implemented climateforced firn densification models. In the following model experiments, we will be using the Herron-Langway model (Herron and Langway, 1980) in its stress-based formulation as the input model. While the HL model on one hand is capable of reproducing 175 the firn densities in the investigated NEGIS region accurately, its stress-based formulation additionally considers potential strain-induced changes of the overburden load, which an accumulation-based formulation does not capture.

Data
Modeling strain softening requires knowledge about the horizontal strain rates. For the modeling experiments conducted in this paper, horizontal strain rates are computed from surface velocity maps of the Greenlandic (GrIS) and Antarctic ice sheet (AIS) using the logarithmic strain rate computation method as discussed by Alley et al. (2018). Nominal computation of the strain rates proved however to be sufficient as difference between the two computation methods were smaller than the uncertainty induced by the velocity data itself. Before determining the strain rates from the velocity fields, a Gaussian filter is applied on the velocity maps to reduce the impact of processing artifacts in the data, which likely were caused by combining velocity data from different sources for producing these velocity fields. These artifacts appear as a grid structure in the strain rate products and clearly do not 190 represent any physical information, but would lead to an overestimation of the horizontal strain rates, if not removed. In the GrIS velocity data a Gaussian filter with a standard deviation of 2 pixels was applied. For the AIS we use a variable Gaussian filter with standard deviations between 2 − 10 pixels. Regions with poor data coverage and higher reported uncertaintiese.g. in the polar hole -are smoothed more. The smoothing reduce the effective spatial resolution of the velocity grid, which can hide the maximum local strain in shear margins. However, this is not a major concern as km-scale variations of horizontal 195 strain rates average out over the firn age.
For validating the strain softening model, firn density data recorded at the NEGIS in the vicinity of the EGRIP ice core site is used. As suggested by Vallelonga et al. (2014) the NEGIS with its high-strain shear margins offers excellent opportunities for studying firn densification processes. The data reproduced in this study comprises a 37 km-long cross section of the NEGIS firn densities recorded with active seismic surveying by Riverman et al. (2019) and the density of the NEGIS firn core (Vallelonga 200 et al., 2014). Their locations are shown in Fig. 1. Additionally, the density of the EGRIP S5 2019 shear margin firn core is modeled with the intent to compare the model with directly measured firn density data from a high-strain area. However, this data is unfortunately not available yet, as the firn core is stored at the EGRIP station and is not accessible, because of COVID-19-related restrictions of field work.
The horizontal strain rates that the firn at these sites has experienced in the past, are computed by step-wise backtracing 205 their position according to the velocity field with a monthly resolution and interpolating the computed horizontal strain rate components to these points at every time step. We force the model with a constant temperature of −29.9 • C. This is the seasonality-corrected mean of the 10 m temperature recorded between June 2019 and January 2021 at the PROMICE weather station at EGRIP (Fausto and van As, 2019). Using a constant temperature input is justified, because we are mainly interested in firn processes occurring below a depth of 10 m, where seasonal and inter-annual variability of temperature is smoothed out by 210 heat conduction. Further, we do not expect a significant spatial variability of temperature over this relatively small study region.
As accumulation rate input we use the values derived by Riverman et al. (2019) from the Accumulation Radar of an Operation IceBridge flightline (Paden et al., 2014(Paden et al., , updated 2018 crossing the NEGIS density profile, see Fig. 1b. We extrapolate the accumulation rate from the flight line to the sites where densification is modelled using nearest point interpolation. The effect of strain softening on firn densification is not only studied on local, but also on ice sheet wide scales. For this 215 purpose the output data of the regional climate model HIRHAM5, forced by the ERA-Interim reanalysis product (Dee et al., 2011), is employed as climatic forcing. For the GrIS the mean of the precipitation and surface temperature output between 1980 and 2014 is used Mottram et al., 2017) and has a spatial resolution of 0.05 • , respectively ∼ 5 km.
For the AIS the mean surface mass balance and 10 m temperature output between 1980 and 2017 are used (Hansen et al., 2021), which have a resolution of 0.11 • , respectively ∼ 12.5 km.

220
Newly formed surface layers are set to a density of 315 kg m −3 in all model runs, following Fausto et al. (2018). This matches well with measured surface densities in the NEGIS region (Schaller et al., 2016). Although this assumption is certainly not valid over the whole GrIS and AIS, the potential bias can be neglected in our ice sheet wide experiments, as the main interest lies in the identification of the relative changes in the predicted firn properties between models that do and do not consider strain softening.

Sensitivity test
In order to understand how strain softening behaves under different dynamic conditions a sensitivity test was conducted. For the climatic conditions present at the EGRIP ice core site in the center of NEGIS the firn density and firn age was modeled for a range of effective horizontal strain rates betweenε E,H = 0 andε E,H = 7 × 10 −3 yr −1 .

230
The strain dependent age profile of the firn is shown in Fig. 2a. The picture is similar for the sensitivity of the firn age at a certain depth and for the BCO age, represented in Fig. 2a by the age of the firn at the BCO density line. Again both quantities are most sensitive to the effect of strain softening at low strain rates. The shift of the firn age at a given depth is however rather weak with a decrease of the firn age by up to 10 % for 240ε E,H = 7 × 10 −3 yr −1 . The BCO age in contrast is strongly affected by the horizontal strain rates with a decrease of more than 50 % at the highest strain rate values in this test.
In summary, both the BCO depth as well as the BCO age are strongly affected by high strain rates, but even low strain rates of less than 1 × 10 −3 yr −1 do affect the firn properties considerably.

Comparison with firn cores 245
The firn density profiles at the sites of the NEGIS firn core and the EGRIP S5 2019 firn core are modeled with the HL model by first considering no strain and subsequently by additionally activating the modules for horizontal divergence (not shown), strain softening and finally adding a further correction for strain softening, which will be introduced in the following.
In Fig. 2b the results for the NEGIS firn core, drilled at the site of the EGRIP station, are shown and compared to the data.
As the effective horizontal strain rate is small at this site, the differences between the various model setups are also small. The 250 no strain-model already reproduces the density data with good agreement. While horizontal divergence only contributes to a reduction of the firn thickness by 1 m, strain softening thins the firn by additional 6 m.
As the no-strain model already matches the data, the strain softening model underestimates the firn thickness. This can be explained by the fact that the no-strain HL model has been empirically tuned to the density profiles at a range of sites. Thus, it implicitly also accounts for some effective strain rate that is typical for these sites. The mean effective horizontal strain 255 rate at the sites of the firn cores, used to calibrate the HL model, amounts toε cor = 4 × 10 −4 yr −1 . Hence, a corresponding strain softening contribution can be expected to be captured by the HL model, which in regard of the sensitivity test results is considerable.
We therefore apply a correction that reduces the effective horizontal strain rate input into the strain softening model, by subtracting a value ofε cor from the effective horizontal strain rate, respectively by setting lower strain rates to 0, as the 260 effective horizontal strain rate may not become negative. Becauseε cor almost matches the effective horizontal strain rate at the EGRIP site, the correction nearly completely rules out the effect of strain softening for this site and only the contribution of horizontal divergence is what remains.
Although the shear margin firn density data of the S5 2019 firn core is not available, the corresponding modeled firn density profiles are shown in Fig. 2c as an example for a high-strain site. The modeled density profile for the cases of no strain 265 resembles the profiles at the EGRIP site. While horizontal divergence again only slightly alters the firn profile, strain softening causes a significant reduction of the firn thickness by 28 m and thereby causes the kink at the critical density to disappear. At this high-strain site the previously introduced correction factor only has little affect on the firn density profile, which is caused by the lower sensitivity of the firn densification to changes of the effective horizontal strain rate at high values.

Validation with NEGIS firn density cross section 270
As a next step the firn densities recorded along a cross section of NEGIS by (Riverman et al., 2019) are reproduced. The original data is shown in Fig. 3a with the firn-ice transition being indicated by the white contour line. In this profile the increased firn density within the shear margins of NEGIS can clearly be seen, causing a reduction of the firn thickness by 20 to 30 m.
In Fig. 3b the density profile is modeled with the corrected strain softening model at the locations of the individual seismic survey sites as shown in Fig. 1b, according to the forcing parameters displayed in Fig. 3c, which are the accumulation rateḃ 275 and the effective horizontal strain rates experienced along the flow path, illustrated as the mean effective horizontal strain rate over the firn age,ε E,H . For the temperature a constant value of −29.9 • C was assumed. Additional to the results of the corrected strain softening model and its BCO contour, the latter is also displayed for the cases of no strain, horizontal divergence, and the uncorrected strain softening model.
In the case of no strain, the BCO line represents the impact of the accumulation rate variability alone. It shows that the ob-280 served density peaks in the shear margins cannot be attributed to the accumulation pattern. In contradistinction to the observed lower firn thickness, the higher snow accumulation in the shear margin would even promote an increase of firn thickness.
Similar conclusions need to be drawn for the effect of horizontal divergence, as the corresponding BCO line only differs slightly from the case of no strain. Horizontal divergence merely affects the firn thickness by a few meters, whereat no clear trend for a firn thinning can be seen, but instead firn thickness also increases where velocities are actually converging as it is 285 for example the case at the S5 2019 firn core site, indicated by the left vertical line. If velocities diverge at one place, they must converge elsewhere. Therefore, horizontal divergence cannot explain a pure firn thinning pattern, but always results in an increase of firn thickness nearby and consequently it particularly cannot explain the reduced firn thickness in the shear margins of ice streams.
The increased firn density in the shear margins with the respective lowering of the BCO depth can only be reproduced, when 290 strain softening is included into the model. In this case the extent of the density peaks can be reproduced well, which validates the model and supports the idea that the enhanced firn densification rates at high strain rates are caused by strain softening. Comparing the strain softening cases with and without subtracting the correction factor ofε cor = 4 × 10 −4 yr −1 along the cross section, it is becomes notable, that even such a small contribution can considerably alter the modeled firn thickness by up to 5 m. As the sensitivity of firn densification to strain softening is highest at low strain rates, the correction also has the 295 biggest impact in-between and outside of the shear margins, where strain rates are low. Within the high-strain shear margins the correction on the contrary has only very little effect.
The main difference between data and model lies in an apparent shift of the modeled firn density profile of 2 km towards south-east. No, processing error that would explain such a shift could be identified, neither for the model nor the data. Instead, the shift could point towards a potential movement of the ice stream position over time, which would not be represented in 300 the model as the past horizontal strain rates are inferred from present day surface velocities, whereby the assumption of steady state conditions is implicitly made.

Firn properties along NEGIS
In Fig. 4a the modeled change of the firn thickness is compared with the elevation of the ice sheet surface along the seismic survey line. The shear margin troughs resemble the modeled reduction of the firn thickness with regard to depth, extent and 305 location, suggesting that they are formed by a collapse of the firn layer due to strain softening. However, as the agreement is not perfect, other factors, such as upstream effects, accumulation variations and the sub-glacial topography, must also alter the surface elevation and the structure of the shear margins troughs.
Previously, the expected depth of the shear margin troughs was estimated by integrating the vertical strain rate caused by horizontal divergence along the flow line, giving a total strain of −0.1 which translates into a trough depth of 200 to 300 m 310 (Fahnestock et al., 2001). Our results in contrary indicate that horizontal divergence in the firn does not contribute to the trough formation. Although our study focuses on the firn layer, this conclusion can be extended to the ice layer below. Because the firn thinning matches the depression of the surface topography, the ice layer itself has an approximately flat surface, suggesting that no thinning is occurring within. Instead, it indicates that the pattern of horizontal strain rates must be changing with depth. This depth variable horizontal strain rate in turn might have an impact on strain softening. The presented model however assumes 315 that horizontal strain rate are constant with depth, the role of its depth variability hence needs to be investigated in future studies.
The BCO age along the profile is shown in Fig. 4b. While the firn reaches an age of up to 400 yr in the center of the ice stream, the firn-ice transition in the shear margin already takes place after about 200 yr. This means that in the shear margins the densification rate is effectively doubled by strain softening. The BCO age is reduced by 50 % over a distance of only 5 km, 320 whereas the climatic forcing on firn air processes can be expected to vary very little over such short distances. For this reason, we suggest to exploit the shear margins as a natural laboratory for firn air studies, because Characteristic parameters of firn air processes can potentially be better constrained by exploiting this setup.

Implications for ice core dating
The gas enclosed in bubbles at the BCO is younger than the surrounding ice. This introduces a difference between the age 325 of the ice matrix and the gas in an ice core (∆age). Accurate models of densification are thus important for synchronizing ice core records. A precise relative timing is necessary to understand the sequence of events and the physical mechanisms behind past changes in climate (Pedro et al., 2012). E.g. Buizert et al. (2015) finds that abrupt Greenland warming events lead corresponding Antarctic cooling onsets by 218 ± 92years. This conclusion hinges on the estimated ∆age at the WAIS Divide ice core. In contrast to this, Svensson et al. (2020) finds a 122 ± 24year lag between Greenland and Antarctic ice core 330 records using a volcanic synchronization that does not rely on densification processes. To gauge the potential impact of strain softening at this site we test the sensitivity to an effective strain rate of 1 × 10 −3 yr −1 using WAIS Divide climate conditions.
For Holocene climate we find that strain softening reduces the BCO age by 23 % (from 309 to 238 years) and BCO depth by 20 %. For Last Glacial Maximum conditions (−41 • C and 0.1myr −1 ) the BCO age is reduced by 33% or 209 years and BCO depth by 29 %. The impact of strain softening therefore depends on the climatic forcing at the site and can alter over time even 335 if the effective horizontal strain rate remains constant. For the WAIS Divide site the BCO age reduction by a moderate strain rate decreased from 209 to 71 years between the Last Glacial Maximum and the Holocene. This decrease is therefore on the order of the observed time lag between Greenland and Antarctic ice core records and needs to be considered for synchronizing ice core records by gas isotopes.
In classical climate-forced densification models the densification rate and thus ∆age is almost entirely determined by sur-340 face temperature and accumulation rate. Buizert et al. (2021) exploit this to infer past surface temperature from estimates of ∆age and accumulation rate. However, our modelling shows that horizontal strain from large scale ice flow lead to enhanced densification rates and should also be taken into account. A complication when modelling past densification rates is that large scale ice flow could have changed over time.

Firn properties of the polar ice sheets 345
Finally, the ice sheet wide impact of strain softening on firn densification is studied. For this purpose, steady state firn density profiles were modeled using the HL model with the strain softening extension for various combinations of forcing parameters that cover the range of climatic conditions and effective horizontal strain rates that are present over the GrIS and the AIS according to the HIRHAM5 output data and the satellite based velocity fields products.
For Greenland, temperature was altered between −29 • C and −17 • C in steps of 2 • C, accumulation rate was logarithmically 350 increased in 7 steps from 75 mm/yr to 1 m/yr and the effective horizontal strain rate was increased in steps of 1 × 10 −3 yr −1 from 0 to 7 × 10 −3 yr −1 , giving 392 different combinations of forcing parameters in total. For Antarctica temperature was linearly increased from −60 • C to −10 • C in 6 steps, accumulation rate was again logarithmically increased in 9 steps from 5 mm/yr to 1 m/yr and the effective horizontal strain rate was increased up to 10×10 −3 yr −1 with the same spacing as before.
Which gives 594 different combinations for the AIS.
Locations with warmer temperatures and lower accumulation rates than given by these input ranges where not modeled and for the GrIS also places with an average annual melt of more than 1 mm were excluded. This was done, because in the ablation zone of the polar ice sheets additional processes like melt-water percolation and refreezing contribute to firn densification, which will not be enhanced by strain softening. Drawing conclusions about how strain softening affects the general firn densification process in these areas is therefore not easily possible. To not overestimate the contribution of strain 360 softening in these areas, we decided to restrict our studies to the dry zone of the polar ice sheets by introducing the above restrictions.
The BCO depth and BCO age for the input forcing where then determined from the modeled steady state firn profiles. Local firn properties at every point on the ice sheet were obtained by linear interpolation of the local climatic forcing to the parameter grid. With this approach the ice sheet wide contribution of strain softening to firn densification can be studied by comparing 365 the interpolated firn thickness in the cases of no strain to the case when the uncorrected strain softening model is employed.  Fig. 5f, it can moreover be noted that even in the interior of the East Antarctic Ice Sheet (EAIS) strain softening enhances firn densification by up to 10 %, despite of low flow velocities and therefore also low horizontal strain rates being present. This unexpected observation can be explained by very low temperatures and accumulation rates occurring in this region, which give rise to extremely low firn densification rates with a BCO age on the order of 10 3 yr and in consequence enable strain softening to have a relatively strong impact by being active over a long time span.

375
Strain softening hence contributes to firn densification over wide areas of the ice sheets. This however does not mean that existing firn densification models, which did not consider strain softening, generally overestimated firn thicknesses, but rather that some contribution of strain softening might have been captured by the model parameters representing the effect of temperature and accumulation rate. The sensitivity of those classical models to the climatic forcing will hence not be accurate.
To illustrate this effect, Fig. 5 additionally shows the locations of the firn cores that were used for empirically tuning the  Langway model. At least at some of these sites, especially on the AIS, firn thickness is appreciably affected by strain softening.
The mean effective horizontal strain rate over all sites where data is available amounts to 4 × 10 −4 yr −1 . Here, the site of the Little America V firn core, which Herron and Langway (1980) already noted to be affected by horizontal stress, is not even included. Consequently, it can be expected that the empirically tuned model parameters for the temperature and accumulation rate dependence implicitly considered strain softening corresponding to an effective horizontal strain rate of 4 × 10 −4 yr −1 .

385
To take account of this implicit contribution, the previously introduced correction factor was included into the model. However, the correction only works properly at effective horizontal strain rates aboveε cor . Below this values the input strain rates can only be reduced to zero, because the effective horizontal strain rate must not be negative. Hence, the correction is not suitable for compensating for the strain softening effect in the HL model in areas with very low strain. To accurately represent to consider them together during the empirical tuning of the firn model. We leave this to future studies.
Strain rates derived from remotely sensed velocities are sensitive to the degree of smoothing applied. Spatially uncorrelated velocity noise will lead to a positive bias in the effective strain rate. Smoothing reduce the noise amplitude, and will act to lessen this bias. Unfortunately, smoothing will also blur the true strain rate signal leading to a negative bias in the effective strain rate in shear margins. The degree of smoothing is therefore a compromise between reducing noise, and not degrading the 395 signal too much. In this paper, we have chosen the degree of smoothing necessary to remove obvious artefacts. However, in the interior regions with very slow flow the true strain rates may be so small that even a tiny remaining noise amplitude can still be a substantial component of the estimated effective strain rate. This is an important caveat when interpreting the continental scale maps in fig. 5.

400
We have developed a extension for firn densification models that is capable of correcting the densification rate of any climateforced firn model for the effect of strain softening. Employing this model, it was studied how strain softening affects firn densification on local and ice sheet wide scales.
We found that the sensitivity of firn densification to strain softening is highest at low strain rates and that therefore even low strain rates can affect the firn thickness considerably in areas where forcing by accumulation rate and temperature is weak. In 405 high-strain areas, such as the shear margins of ice streams, a significant acceleration of firn densification by strain softening was modeled, which is in good agreement with observed lower firn thickness in these areas. As other potential processes, like horizontal divergence or a greater accumulation in the shear margin troughs, could not explain this reduction of firn thickness, our work supports the idea that strain softening is the principle cause. It was further observed that the change of firn thickness resembles the lowering of the surface elevation in the shear margins, which suggests that the shear margin troughs form because 410 of a faster settling of the firn due to strain softening.
Strain softening not only affects firn thickness, but also reduces the age of the firn at the firn-ice transition. According to our model this can lead to a reduction of the BCO age by around 50 % over a few kilometers in the shear margin of ice streams.
We therefore suggest to exploit this feature as a natural laboratory in future firn air studies, because climatic forcing over such small distances will only vary slightly. Moreover, strain softening can strongly alter the BCO age over time, even at constant, 415 moderate strain rates. For ice core dating this induces a bias in the BCO age, which previously has not been considered, but is on an considerable order for synchronizing ice core records by gas isotopes.
Finally, we demonstrate that strain softening has a substantial effect to firn densification over wide areas of ice sheets and as a consequence that horizontal strain rates should generally be considered in firn densification modeling, because a restriction to climatic forcing parameters results in a misrepresentation of the latter. Our work therefore suggests that besides 420 temperature and accumulation rate also the effective horizontal strain rate should be considered as a relevant forcing parameter