Firn densification in two dimensions: modeling the collapse of snow caves and enhanced densification in ice-stream shear margins

Journal of Glaciology Article Cite this article: Arrizabalaga-Iriarte J, Lejonagoitia-Garmendia L, Hvidberg CS, Grinsted A, Rathmann NM (2025) Firn densification in two dimensions: modeling the collapse of snow caves and enhanced densification in ice-stream shear margins. Journal of Glaciology 71, e61, 1–16. https:// doi.org/10.1017/jog.2025.6 Received: 4 April 2024 Revised: 6 January 2025 Accepted: 13 January 2025 Keywords: Arctic glaciology; Ice physics; Ice rheology; Ice streams; Polar firn Corresponding author: Jon Arrizabalaga-Iriarte; Email: [email protected] © The Author(s), 2025. Published by Cambridge University Press on behalf of International Glaciological Society. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/ by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited. cambridge.org/jog Firn densification in two dimensions: modeling the collapse of snow caves and enhanced densification in ice-stream shear margins Jon Arrizabalaga-Iriarte1,2,3, Lide Lejonagoitia-Garmendia2, Christine S. Hvidberg2, Aslak Grinsted2and Nicholas M. Rathmann2 1Basque Centre for Climate Change (BC3), Leioa, Spain; 2Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark and 3Faculty of Science and Technology, University of the Basque Country (UPV/EHU), Leioa, Spain Abstract Accurate modeling of firn densification is necessary for ice core interpretation and assessing the mass balance of glaciers and ice sheets. In this paper, we revisit the nonlinear-viscous firn rheology introduced by Gagliardini and Meyssonnier (1997) that allows multidimensional firn densification problems to be posed, subject to arbitrary stress and temperature fields. First, we extend the calibration of the coefficient functions that control firn compressibility and viscosity to five additional Greenlandic sites, showing that the original calibration is not universally valid. Next, we demonstrate that the transient collapse of a Greenlandic firn tunnel can be reproduced in a cross-section model, but that anomalous warm summer temperatures during 2012–14 reduce confidence in attempts to independently validate the rheology. Finally, we show that the rheology can explain the increased densification rate and varying bubble close-off depth observed across the shear margins of the Northeast Greenland Ice Stream. Although we suggest more work is needed to constrain the near-surface compressibility and viscosity functions of the rheology, our results strengthen the empirical grounding of the rheology for future use, such as modeling horizontal firn density variations over ice sheets for mass-loss estimates or estimating ice-gas age differences in ice cores subject to complex strain histories. 1. Introduction Snow precipitated over glaciers and ice sheets transforms into ice in the uppermost layers of the ice column through firn densification. During this process, old snow is buried and compressed as a result of the overburden pressure until it reaches the density of pure ice. Understanding and modeling the details of this process is central for the interpretation of ice core records and assessing the mass balance of glaciers and ice sheets based on satellite altimetry. In the case of ice core records, a good understanding of the air-trapping mechanism is key for interpreting gas records. The air present in firn pores becomes isolated from the atmosphere not until deep in the firn column, at the so-called bubble close-off (BCO) depth. This can create 100to 1000-year differences in age between trapped gases and the surrounding ice (Δage), which must be considered to properly synchronize records (Schwander and others, 1997). In the case of mass-balance monitoring from satellite altimetry, the problem is to translate changes in altitude into changes in mass. Sea-level rise estimates can be biased if the firn density field (i.e. firn compaction) is not correctly accounted for (Sørensen and others, 2011, Lipovsky, 2022), and recent efforts to include horizontal density variations over the West Antarctic ice sheet in large-scale modeling finds a correction in volume above flotation over 40 years of 10% (Schelpe and Gudmundsson, 2023). Firn densification models are an important tool that can help estimate the effect of densification on ice-core climatic records and mass-loss uncertainties. Ideally, densification models would be derived from first principles, but densification is a complex process in which crystals rearrange and change their shape and size. As a consequence, the relative importance of the different mechanisms driving densification is divided into several stages depending on local conditions (for details see Cuffey and Paterson, 2010, p. 16–22). To simplify matters, large-scale models of firn columns therefore tend to be constructed based on phenomenological arguments. Herron and Langway (1980) proposed a one-dimensional empirically tuned model of the first and second stages of densification for predicting depth–density, depth–load and depth–age profiles based exclusively on local temperature and accumulation rate. Their celebrated model has since become the benchmark for comparing more sophisticated models with (e.g. Barnola and others, 1991, Li and Zwally, https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press 2 Jon Arrizabalaga-Iriarte et al. 2011). However, reproducing observed density profiles at sites with diverse climatic conditions remains challenging, and some models give inconsistent results for the same boundary conditions (Lundin and others, 2017), suggesting that more work is warranted. Densification models that attempt to expand on the physical rigor of previous work have also been proposed, although improved accuracy over empirically based models does not necessarily follow (Thompson-Munson and others, 2023). For example, Alley (1987) proposed a simple model for porous firn densification by grain boundary sliding, several models make use of a one-dimensional compactive viscosity (Arnaud and others, 2000, Morris and Wingham, 2014, Stevens and others, 2023), and threedimensional constitutive relations exist that argue for certain viscosity functions of firn (Salamatin and others, 2009, Fourteau and others, 2024). Finally, we note that Arthern and others (2010) proposed a Nabarro–Herring creep equation for porous materials that is the basis for several of the models used for altimetry corrections (Smith and others, 2020), and SNOWPACK—a physically based land-surface snow model originally developed to support avalanche warning (Bartelt and Lehning, 2002, Lehning and others, 2002a,2002b)—has also been successfully adapted to model near-surface firn densification (Keenan and others, 2021). 1.1. Beyond one-dimensional modeling A central problem is to take into account the larger-scale ice flow, as there is evidence of accelerated firn densification when horizontal strain rates are nonzero (Kirchner and others, 1979, Alley and Bentley, 1988, Riverman and others, 2019). This phenomenon, called strain softening, was first introduced by Alley and Bentley (1988); they argued that during the second stage, densification is primarily driven by dislocation creep, which increases with the square of the effective stress. Hence, in regions where nonnegligible horizontal strain rates are superimposed on the vertical firn compaction, the effective viscosity is reduced, which should, in turn, enhance densification. Motivated by this, Oraschewski and Grinsted (2022) recently extended the Herron and Langway (1980) model (HL) by introducing a scale factor that allows the inclusion of horizontal strain rates as a sort of forcing parameter, together with the commonly used climate variables (temperature and accumulation rate). Specifically, the densification rate derived from the classical HL model is multiplied by a scale factor that is a function of the otherwise unresolved strain rate components. Although extensions based on Herron and Langway (1980) are generally computationally inexpensive, they remain one-dimensional. Multidimensional approximations can be constructed by horizontally joining firn-column models (e.g. Oraschewski and Grinsted, 2022), but formulating firn densification in terms of a three-dimensional constitutive stress–strain rate relationship is desirable for several reasons: (i) the effect of nonzero horizontal strain rates can be represented naturally (without the need for any ad hoc inclusion), (ii) twoor three-dimensional firn compaction problems can be more easily solved for complicated boundary conditions and geometries using, e.g. the finite element method and (iii) Glen’s flow law for solid ice may be extended to also characterize porous firn, thereby allowing for a seamless numerical treatment of the entire firn–ice column. The semi-empirical rheology proposed by Gagliardini and Meyssonnier (1997)—henceforth GM97—is exactly such a field formulation. It is based on a general compressible power-law rheology for porous materials, adapted to fit in situ measured density profiles from Site-2, Greenland (see Fig. 2), and cold room deformation tests. The GM97 rheology has previously been shown capable of reproducing age–depth profiles of mountain glaciers (Lüthi and Funk, 2000, Zwinger and others, 2007, Gilbert and others, 2014, Licciulli and others, 2020) and predicting the evolution of snow caves buried over time in Antarctica (Brondex and others, 2020). The GM97 rheology depends on two coefficient functions of density that determine the local compressibility and viscosity. The functional form of these, and associated calibration constants, have received little attention in the literature apart from the initial treatment by Gagliardini and Meyssonnier (1997) and subsequently by Zwinger and others (2007). The proposed coefficient functions mainly disagree for near-surface firn densities, which can have a large impact on modeled surface velocities and densification rates (Gilbert and others, 2014). This has led to some disagreement in the literature about which coefficient functions to use. More work is therefore warranted to make clear whether the calibrations proposed by Gagliardini and Meyssonnier (1997) and Zwinger and others (2007) (i) hold for a wider range of sites, where firn cores have since been drilled and (ii) hold in two-dimensional settings, where nontrivial densification fields have been observed. If so, this would put the GM97 rheology on stronger empirical grounds and increase confidence in its future use, which is the primary aim of our work. In this paper, we test the flexibility and usefulness of the GM97 rheology for modeling firn compaction beyond one-dimensional problems. First, we revisit the original calibration of the rheology and assess the universality of it by trying to reproduce five additional Greenlandic ice core sites, where one-dimensional modeling is appropriate. Then, we test the model’s accuracy and dependency on temperature by attempting to reproduce the collapse of a near-surface snow cave, constructed at the NEEM ice core site in Greenland. Finally, we test whether the rheology can reproduce the enhanced densification observed at the shear margins of the Northeast Greenland Ice Stream (NEGIS), believed to be caused by strain-softening effects. In the following, we begin by introducing the GM97 rheology before considering the special oneand two-dimensional cases thereof, needed for our experiments. 2. Method For a compressible material, mass conservation implies that 𝜕𝜌 𝜕t+∇⋅(𝜌u)=0,(1) where uand 𝜌are the velocity and density fields, respectively. The velocity field of a Stokes fluid is governed by the momentum balance ∇⋅𝝈+𝜌g=0,(2) where 𝝈(. 𝝐)is the stress tensor, . 𝝐=(∇u+(∇u)T)/2 is the strain rate tensor and gis the gravitational acceleration vector. The GM97 rheology treats firn densification as a secondary creep problem involving u,𝜌and the temperature (or enthalpy) T, assuming negligible effects from primary creep, snow metamorphism and the brittle fracturing of snow. The rheology represents both firn and ice using a single, compressible Norton–Bailey power-law fluid that conforms to Glen’s incompressible flow law in the limit 𝜌→𝜌ice, where 𝜌ice =917 kg m−3is the density of glacier ice. The rheology is to be understood as an extension of Glen’s flow law that depends on the first invariant of the strain rate tensor, https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press Journal of Glaciology 3 tr(. 𝝐), in addition to the usual second invariant, . 𝝐:. 𝝐. The rheology was first proposed by Duva and Crow (1994) for general porous materials but has since been adopted in glaciology, and recently with renewed interest (Brondex and others, 2020, Licciulli and others, 2020). Written in its inverse form (see Appendix B) required by the momentum balance (2), the rheology is 𝝈=A−1/n. 𝜖(1−n)/n E[1 a(. 𝝐−tr(. 𝝐) 3I)+3 2btr(. 𝝐)I],(3) . 𝜖2 E=1 2a(. 𝝐:. 𝝐−tr(. 𝝐)2 3)+3 4btr(. 𝝐)2,(4) where n=3 following the above-mentioned adoptions of the rheology, . 𝜖Eis the effective strain rate and aand bare the coefficient functions controlling firn viscosity and compressibility (elaborated on below). Note here that Glen’s law is recovered when a=1, b=0 and tr(. 𝝐)=0, which must be fulfilled in the limit 𝜌→𝜌ice. Since the flow-rate factor, A, depends on temperature, the evolution of the temperature field (T) needs to be considered, too: 𝜌c(𝜕T 𝜕t+u⋅∇T)=∇⋅(kT∇T)+𝝈:. 𝝐, (5) where c=c(T)and kT=kT(T,𝜌)are the specific heat capacity and thermal conductivity of ice (see Appendix A for empirical functional forms used). In GM97, the flow-rate factor is modeled the usual way as an Arrhenius activation function of temperature (Greve and Blatter, 2009, Greve and others, 2014): A(T)=A0exp[−Q/(RT)], (6) where, following Zwinger and others (2007), the canonical dependence on the temperature of the prefactor A0and the activation energy for creep, Q, are presumed: A0={3.985 ×10−13 s−1Pa−3for T≤−10∘C 1.916 ×103s−1Pa−3for T>−10∘C,(7) Q={60 kJ mol−1K−1for T≤−10∘C 139 kJ mol−1K−1for T>−10∘C.(8) The fluidity also depends on other factors such as pressure (Weertman, 1983), water content (Barnes and others, 1971, Duval, 1977), impurities (H orhold and others, 2012) and grain sizes and orientations (Duval, 1973, Shoji and Langway, 1988). The pressure dependence is estimated to be relatively small, even at hydrostatic pressures near ice sheet beds, and therefore neglected (Rigsby, 1958). Regarding water content, temperatures are generally considered to be low enough to neglect the presence of liquid water. Impurities, such as dust, are known to increase ice softness by up to a factor of two, but in the Greenland firn column considered here— which consists only of Holocene deposition—the dust loading is minimal or can at least be assumed to be uniform to first order. Finally, changes in the size and orientation of grains are thought to be important in dynamic regions such as ice streams and their shear margins (e.g. Gerber and others, 2021), but expanding (3)–(4) to also allow for viscous anisotropy is a daunting task and out of scope of this work, so the effect of developed crystal fabrics in the firn constitute a rheological uncertainty throughout this work. Figure 1. Coefficient functions a(red) and b(blue) proposed by Gagliardini and Meyssonnier (1997) and Zwinger and others (2007) for n=3. 2.1. Coefficient functions Following Duva and Crow (1994), the coefficient functions aand bare taken to depend exclusively on the relative density 𝜌= 𝜌 𝜌ice .(9) In this way, if aand bare affected by temperature, grain sizes, etc., it is implicitly assumed that such dependencies can be factored out and absorbed into A. Indeed, this separation of dependencies is desirable since (3)–(4) can then be easily made to reduce to Glen’s law in the limit 𝜌→1, as noted above. The coefficient functions aand bdepend in part on a model for how stresses concentrate in the presence of enclosed voids (air) (Duva and Crow, 1994) for densities above the critical value 𝜌≥ 𝜌crit =0.81: a0=1+2 3(1− 𝜌) 𝜌2n/(n+1),(10) b0=3 4(n−1(1− 𝜌)1/n 1−(1− 𝜌)1/n)2n/(n+1),(11) and in part on an empirical scaling relation for densities 𝜌<𝜌crit that has been found to fit cold room experiments and densification measured at Site 2, Greenland. Specifically, two such sets of aand b scalings have been proposed: Gagliardini and Meyssonnier (1997) suggested (Fig. 1, dashed lines) a1=(a0/b0)b1(12) b1=⎧ { { ⎨ { { ⎩ exp(451.63 𝜌2−474.34 𝜌+128.12)for 𝜌<0.5 exp(−17.15 𝜌+12.42)for 0.5≤ 𝜌≤ 𝜌crit b0for 𝜌crit <𝜌 (13) https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press 4 Jon Arrizabalaga-Iriarte et al. whereas Zwinger and others (2007) suggested (Fig. 1, dark solid lines) a1={kaexp(−𝛾a( 𝜌− 𝜌sfc)) for 𝜌≤ 𝜌crit a0for 𝜌crit <𝜌 (14) b1={kbexp(−𝛾b( 𝜌− 𝜌sfc)) for 𝜌≤ 𝜌crit b0for 𝜌crit <𝜌 .(15) Here 𝜌sfc =0.4 is the assumed relative density of surface snow, continuity at 𝜌= 𝜌crit implies the scaling exponents are 𝛾a=ln(ka/a0( 𝜌crit)) 𝜌crit − 𝜌sfc and 𝛾b=ln(kb/b0( 𝜌crit)) 𝜌crit − 𝜌sfc ,(16) and kaand kbare parameters to be calibrated against observations or experiments, suggested by Zwinger and others (2007) to be k≡ka=kb≃1000.(17) Note that kaand kbare the values of aand bat the surface, 𝜌= 𝜌sfc. As pointed out by Brondex and others (2020), kaand kbshould ideally be re-calibrated on a case-by-case basis, although adopting the cold room and Site 2 calibration proposed by Zwinger and others (2007) (k≃1000) has been found sufficient in several cases (Gilbert and others, 2014, Brondex and others, 2020, Licciulli and others, 2020) and is regarded as the canonical value. In this work, we revisit the best-fit value of kby considering additional oneand two-dimensional model experiments that are evaluated against observations. 2.2. Numerics In the numerical experiments that follow, we solved the coupled density, momentum and thermal problem using FEniCS (Logg and others, 2012), relying on Newton’s method to solve nonlinearities. For reasons explained below, the ice stream scenario is not thermally coupled, but the mechanical problem is solved using the same method. The Jacobian of the residual forms (required for Newton iterations) were calculated using the unified form language (Alnæs and others, 2015), used by FEniCS to specify weak forms of partial differential equations (PDEs), which supports automatic symbolic differentiation. All weak forms are presented in Appendix A. For our two-dimensional experiments, meshes were constructed using gmsh (Geuzaine and Remacle, 2009) and updated between time steps to evolve the interior and exterior free-surface boundary. 3. Validation of rheology We considered three independent numerical experiments to test the proposed best-fit value of the calibration constant, k≃1000 (Zwinger and others, 2007). First, we swept over kto determine which value can best reproduce observed firn density profiles from a wider range of Greenlandic ice core drill sites (case 1). Second, we sought a value of kthat could best reproduce the measured collapse of a trench constructed at NEEM (Steffensen, 2014) (case 2). Third, we attempted to reproduce the observed enhanced densification over the shear margins of NEGIS by varying k(case 3). 3.1. Case 1: Reproducing Greenlandic ice cores We considered a total of six Greenlandic firn cores from the deep drilling sites at DYE-3, GRIP, NGRIP, NEEM, Site 2 and Figure 2. Sites of Greenlandic ice cores used to validate the GM97 rheology and determine k. Satellite-derived velocities from the MEaSUREs program are shown in colored contours (Howat, 2020). Site A (Fig. 2). Since firn temperatures are approximately depthconstant at each site (Bréant and others (2017); Table 1), the problem reduces to that of solving for 𝜌and the vertical velocity uz if the density and velocity fields are approximated as horizontally homogeneous (similar to traditional one-dimensional modeling, although this assumption might be less well justified at dynamic sites). Unlike the transient two-dimensional problems considered in the following, we solved directly for the steady states of (1)–(2) subject to the boundary conditions 𝜌(z=0)=𝜌sfc,(18) uz(z=0)=−. a,(19) uz(z=−H)=−. a𝜌sfc 𝜌ice ,(20) where His the height of the modeled firn column, 𝜌sfc is the surface snow density and . ais the snow accumulation rate at the site (Table 1). The boundary condition (20) assumes that, in steady state, the mass flux into the firn column equals that exiting the https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press Journal of Glaciology 5 Table 1. Depth-averaged temperature, accumulation rate, and surface density of each site considered. Accumulation rates and surface densities follow from Bréant and others (2017). For EGRIP, the accumulation rate and surface density data are retrieved from Karlsson and others (2020) and Schaller and others (2016), respectively. Site Depth average T(∘C) . a (kg m−2a−1)𝜌sfc (kg m−3) Site-2 −25.0*(Langway, 1970) 360 350.1 Site-A (Crête) −29.5(Clausen and others, 1988) 282 321.7 DYE-3 −21.0(Dahl-Jensen and others, 1998) 500 357.0 GRIP −31.7(Johnsen, 2003) 210 367.0 NGRIP −31.5(Dahl-Jensen and others, 2003) 175 299.9 NEEM −28.8(Orsi and others, 2017) 200 307.2 EGRIP −28.0*(Zuhr and others, 2021) 130 290 *Average surface temperature was used due to the lack of a temperature profile. column, where Hmust be taken sufficiently large to guarantee that the bottom of the model domain is pure ice, 𝜌(z=−H)=𝜌ice. We determined the best-fit value of kby sweeping a wide range of values and calculating the resulting model–observation misfit. Specifically, we considered k∈[1;1000]and calculated the rootmean-square-error (RMSE) between the modeled and observed density profiles: RMSE =√ √ √ ⎷1 N N ∑ i=1(𝜌(zi)−𝜌obs(zi))2,(21) where Nis the total number of discrete depth levels at which the misfit is evaluated. Computing the percentage error of each level is also an option (weighing the mismatch with respect to the reference value it is deviating from), but the results obtained in this way are virtually identical (data not shown). Fig. 3 shows the RMSE misfits calculated for the lower-density section of the firn column ( 𝜌<0.8), where kis expected to have the largest influence (see the above section on the coefficient functions). The results suggest a misfit minimum exists somewhere between k=100 and k=500 but with disagreements between the sites; that is, a global best-fit kmight not be strictly applicable across all sites (disregarding effects from seasonal temperature variations not treated here for simplicity). For reference, Fig. 4 shows the model performance in the case of Site 2 for different values of k. Note that the sensitivity to the magnitude of kdecreases for relative densities above 0.8. The underestimation of density at depth is present in all the sites modeled. 3.2. Case 2: Trench collapse at NEEM In an attempt to further characterize the model dependence on k, we searched for alternative validation experiments that involve near-surface (low-density) firn, since the choice of kaffects compressibility the most there. We found that the observed collapse of a cylindrical subsurface tunnel, constructed at the NEEM ice core site, is well-suited for this purpose. The tunnel was constructed to determine whether subsurface trenches, built following the Camp Century snow-blowing casting technique (Clark, 1965), but using large inflatable balloons instead, could serve as drilling and science trenches for deep ice core projects. To judge its feasibility, the tunnel collapse rate was tracked for 3 years after its construction. The initial and final geometries are shown in Fig. 5 for reference, although the pictures are from the northern end of the trench, while we modeled the southern end where measurements were made most consistently. Figure 3. Model–observation RMSE of Greenlandic density profiles for k∈[1,1000]. All the curves keep increasing monotonically from k=1000 and onwards (not shown). Figure 4. Modeled Site 2 density profiles for various k(colored lines) compared to observations (Bréant and others (2017); gray dots). The effect of a given kis most evident near the surface: the higher the value of k, the faster the near-surface densification. For 𝜌>0.8, the model generally underestimates densities at a given depth for all k. The tunnel was constructed by blowing snow out of a large trench and installing an inflatable balloon in the open cavity. After inflating it, snow was blown back into the trench again, burying the balloon in a more closely packed (denser) firn with a typical density of 𝜌trench =550 kg m−3(Brondex and others, 2020, Steffensen, 2022). Finally, after a couple of days, the balloon was deflated and removed, preserving the casted, circular structure of the tunnel. We simplify the trench collapse problem by considering a vertical cross-section model, thus implicitly assuming translational invariance along the tunnel (strictly speaking, only relevant https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press 6 Jon Arrizabalaga-Iriarte et al. Figure 5. NEEM balloon trench geometry 2 months after construction on 7 August 2012 (a) and 3 years later on 27 May 2015 (b). No picture of the trench was available for the model target year, 2014, so 2015 is shown instead. Moreover, pictures show the northern end, whereas we used the trench geometry of the southern end that was most consistently measured during the experiment (no pictures available). Pictures were kindly provided by J. P. Steffensen and reprinted with permission according to the NEEM ice-core project media waiver. Figure 6. Zoom-in of the initial geometry and density field of the NEEM trench model. High-density snow was backfilled into the balloon trench, creating a hardened shell surrounding the tunnel compared to the background density field. for very long tunnels). The initial geometry and density field (Fig. 6) was constructed following the technical report by Steffensen (2014) and subsequently allowed to evolve thermomechanically throughout the duration of the test, approximately 2 years long. Note that the roof includes the reported extra 1 m of backfilled snow, but that the presence of other structures in the trench (connecting tunnels, cabins, etc.) is not considered here, which might change results slightly and are a source of uncertainty. The initial background density field is taken to be equal to the smoothed density profile of the NEEM core (Bréant and others, 2017). The snow accumulation rate is set to the annual average over the trench, which is double the reference measurement at NEEM due to wind-driven excess accumulation (Steffensen, 2014). The initial temperature field is taken to be uniform and equal to the firn column average measured at NEEM, ⟨TNEEM⟩(Table 1). However, during the 5 days that it took to construct the tunnel, both the trench and the snow to be backfilled were exposed to higher-thanusual temperatures (between −3 and −6∘C; Steffensen (2014)). Thus, the initial temperature of the firn that surrounded the tunnel was possibly up to 25∘C warmer than the firn-average-background temperature (most likely less, though). In order to assess the potential underestimation of the initial collapse rate, caused by prescribing too cold firn, we therefore also consider a second scenario in which all the backfilled snow starts out with a temperature of −5∘C. This is just an approximation and does not intend to accurately represent the initial temperature field but rather allows us to gauge the impact that a much hotter trench might have on our results. The model domain has a width of L=20 m and a timeevolving height profile h(x,t)relative to a fixed bottom boundary located at a depth of H=30 m below the initial surface. Both Land Hare large enough to avoid boundary conditions affecting the solution (judged by trial and error). The bottom boundary is fixed by a free slip condition and kept at the site’s depthaveraged firn temperature. At the surface, we thermally force the model by setting the firn surface temperature equal to measurements from the local Greenland Climatic Network weather station, TGCnet(t)(Vandecrux and others, 2023). To estimate the impact of not taking the surface temperature evolution into account, we also ran the model by forcing it with a constant firn surface temperature, set equal to the average measured surface temperature over the modeled period, ⟨TGCnet(t)⟩. A periodic boundary condition is imposed on the left and right sides to close the thermal problem. In summary, the boundary conditions are uz(x,z=−H)=0,(22) T(x,z=−H)=Ticecore(z=−H), (23) https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press Journal of Glaciology 7 Figure 7. Modeled evolution of the NEEM tunnel from 10 June 2012 to 27 May 2014. (a) Evolution of the tunnel height for different values of k(line colors) and thermal scenarios (line styles). (b) Corresponding tunnel cross-sections at the end of the simulation. The larger and smaller black curves in panel brepresent the initial and final measured tunnel dimensions, respectively. All simulations are thermodynamically coupled but differ in the surface temperature boundary condition. Experiments denoted by solid and dash-dotted lines use measured time-evolving surface temperatures, whereas dashed lines denote experiments where the average surface temperature was imposed (resulting in practically isothermal conditions). The hot trench scenario includes an initially hotter-than-average backfilled trench (−5∘C), which aims to be more representative of the real initial conditions. The three grey dots in panel ashow measured tunnel dimensions (Steffensen, 2014). The background red line in panel ashows the smoothed hourly temperature record from the site’s weather station that is part of GC-Net (Vandecrux and others, 2023). The original record contains positive temperatures in the first summer (data not shown due to smoothing). T(x,z=h,t)={TGCnet(t)if variable climate, ⟨TGCnet(t)⟩=−27.7∘Cif average climate. (24) Unlike the steady-state model above, the top surface boundary is allowed to evolve freely, as is the interior tunnel boundary. Both boundaries were updated using the Arbitrary Lagrangian-Eulerian method provided by FEniCS, allowing mesh boundary vertices to be displaced according to the local product between velocity and time-step size. In this method, surface accumulation is taken into account by adding a constant vertical velocity component, . a, to the surface boundary. Note that, in effect, surface accumulation gradually leads to an increase in the average height profile since no mass exits the model domain; the domain bottom will, therefore, also gradually tend towards pure ice. Since we are only interested in the relative deformation of the trench geometry—and not the change in centerof-trench position—this is of no concern and does not affect our results. Due to large density gradients causing numerical instabilities, we assume that the surface accumulation has a density of 𝜌trench rather than 𝜌sfc. The accumulation rate was therefore scaled by a factor of 𝜌sfc/𝜌trench =0.56 to ensure that the total mass deposited (and hence overburden load) is correct. Although this assumption reduces the densification rate of the fresh snow layers, this does not have any direct effect on the mechanical evolution of the tunnel below. However, replacing the snow with a thinner layer of denser (and, thus, less airy) firn may cause the tunnel to be more sensitive to surface temperature variations. We solve the transient problem for the four values k= {100,500,1000,2000}using an Euler time-stepping schemewith a time step of Δt=0.75 days. The results are summarized in Fig. 7, where the tunnel height evolution and the final cross-section are shown for each case. Once again, the smaller kis, the stiffer the firn and, thus, the less the tunnel deforms. For values lower than k=2000, the predicted deformation is too small for all the modeled scenarios. For higher values, the tunnel deforms too much, and the ceiling starts to cave inward (not shown here). The hotter the surrounding firn is, the faster the tunnel is found to shrink, as anticipated (dash-dotted compared to solid lines). There is a delayed response in the collapse rates modeled to variations in surface temperature, which is a consequence of the time required for the thermal perturbation to reach the depth of the tunnel. Additionally, we also note that the sensitivity of the tunnel to surface temperature variations decreases with time as it becomes more isolated below the accumulating snow. 3.3. Case 3: NEGIS shear margin troughs The NEGIS is a coherent structure of fast ice flow, initiating less than 150 km from the ice divide and extending more than 600 km to the coast (Grinsted and others, 2022). Elevation maps obtained from ArcticDEM (Porter and others, 2018) show that, compared to average surface height variations, there are relatively pronounced depressions in the shear margins of NEGIS (Fig. 8), and recent work by Hvidberg and others (2020) has verified the depths of the troughs (depressions) in ArcticDEM using GPS stakes. The origin of the shear margin troughs was revealed in a recent seismic survey by Riverman and others (2019) (survey transect is shown in Fig. 9 as a white line, and the surface height in Fig. 10b as a blue line), finding that there is enhanced densification at the shear margins (blue contours in Fig. 10c). Oraschewski and Grinsted (2022) later showed that these depressions might partly be explained by strain softening, as the horizontal velocity gradients are large there. https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press 8 Jon Arrizabalaga-Iriarte et al. Figure 8. ArcticDEM surface topography upstream of NEGIS, hillshaded using the 3D visualization software Blender. The EGRIP drill site is located at the black ball. Note that North has been rotated 135∘clockwise to make the view of the ice-stream margin most clear, hence flow is toward the bottom part of the figure. Figure 9. Satellite-derived surface velocities around EGRIP from the MEaSUREs program (Howat, 2020) and transect (white line) of the 79 geophones used for the seismic survey performed by Riverman and others (2019). Reproducing the shear-margin troughs in a model of Riverman’s transect therefore makes for a good test case of the GM97 rheology (which accounts for strain softening) and potentially allows for another independent way to estimate the best value of k. In addition to the observed shear-margin troughs (blue line in Fig. 10b), we also include the observed BCO depth (i.e. depth at which 𝜌=830 kg m−3) as a calibration target (white line in Fig. 10c). Note that the ability to model BCO depths might be less sensitive to the value of kand more sensitive to the functional forms of aand b(densities are large at the BCO). We model the transect using a setup that combines the approaches from both the one-dimensional firn column and NEEM trench models. The yz domain considered is L=47.5 km wide and H=200 m tall, including a 5 km buffer on both lateral sides to prevent the boundaries from affecting the interior solution (extension not shown in the figures). The boundary conditions are free-slip on the left and right boundary: uy(y=0,z)=0,(25) uy(y=L,z)=0,(26) whereas a non-zero mass flux is imposed on the bottom boundary that balances the surface-integrated mass flux: uz(y,z=−H)=−. a𝜌sfc 𝜌ice .(27) Here, His taken large enough to ensure that the density at the bottom boundary is that of pure ice. The density field is subject only to the surface boundary condition 𝜌(y,z=0)=𝜌sfc.(28) The surface density is considered to be 𝜌sfc =343 kg m−3 (average of the top 2 m of firn, Schaller and others (2016)) while the accumulation rate is taken to be uniform across the transect for simplicity, . a=130 kg m−2 yr−1 (Karlsson and others, 2020), despite it may be up to a 20%higher in the shear margins due to drift snow trapped by the troughs (Riverman and others, 2019). The surface height profile h(y) was updated following the usual kinematic equation for free surface evolution: 𝜕h 𝜕t=. a+u(s) z−u(s) y𝜕h 𝜕y,(29) where . ais the surface accumulation rate, and u(s) yand u(s) zare the surface velocity components in the yz model domain. When remeshing after updating the surface boundary, the density field was linearly interpolated onto the new mesh. In the absence of any horizontal variation in boundary conditions and initial state, the density field would remain horizontally homogeneous in time (the column would evolve solely due to gravitational compaction). However, additional out-of-plane strain rate components play an important role in the NEGIS transect: the shear margins experience large horizontal x–yshear (red line in Fig. 10a) that can cause significant strain softening, and the trunk might experience a slight extending flow (if any) in the along-flow x-direction (blue line in Fig. 10a). We therefore superimpose the observed . 𝜖xy(y)profile on our y–zmodel domain by assuming it to be depth invariant, while neglecting . 𝜖xx and . 𝜖yy as they are comparatively small. All strain rates were calculated based on satellite-derived velocities from the MEaSUREs program (Howat, 2020) and slightly smoothed to avoid numerical instabilities. Finally, lacking information about the . 𝜖xz shear component, we set it to zero following the shallow shelf approximation, commonly used to model ice stream flow. The effect of the out-of-model-plane component . 𝜖xy on strain softening is included by extending the strain rate invariants that https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press Journal of Glaciology 9 Figure 10. (a) Absolute value of the smoothed horizontal strain-rate components along the NEGIS seismic transect (solid lines), and the effective horizontal strain rate . 𝜖eff by Oraschewski and Grinsted (2022), which includes the upstream history, too. (b) Modeled (red and yellow lines) and GPS-measured (blue line) surface elevation anomaly profile along the seismic transect (Riverman and others, 2019). Solid and dashed lines represent isothermal and hot shear-margin experiments, respectively. (c) Modeled BCO depth profiles (lines) plotted on top of the observed density field by Riverman and others (2019). The white line shows the observed 𝜌=830 kg m−3BCO depth contour, and the violet line shows the BCO depth modeled by Oraschewski and Grinsted (OG22). Vertical solid lines show the modeled shear-margin center points (deepest troughs), and dashed lines show the horizontal extent of the imposed 6∘C temperature anomaly in the shear margins. The along flow dimension, x, is pointing out of the plane. enter the effective viscosity, . 𝜖E, according to (to be consistent with) their three-dimensional definitions: tr(. 𝝐)=tr(. 𝝐2D), (30) . 𝝐:. 𝝐=. 𝝐2D :. 𝝐2D +2. 𝜖2 xy,(31) where . 𝝐2D is the two-dimensional strain rate tensor in our yz model. For the temperature field across the transect, we consider two different scenarios. In the first scenario, we assume a uniform temperature of T= −28∘C (and thus a uniform flow-rate factor) (Table 1). In the second scenario, we impose a 6∘C elevated temperature in the ice-stream shear margins as reported by Holschuh and others (2019) (high-end estimate, could be lower), postulated to be caused by strain heating. Since ice is a good thermal insulator, the 6∘C anomaly is specified as an abrupt but depthconstant step function crossing into the shear margins. Although adding a thermal coupling to this problem as well (i.e. evolving the temperature field) would increase the physical realism, the out-ofplane heat flow is poorly known, making the thermal problem not well-posed. On the other hand, by imposing the above-mentioned steady temperature fields, the effect of dynamic strain softening (due to nonlinear viscosity) is made clearer, as its effect can be understood in isolation. For brevity, we report only on four transient simulations that are chosen to be representative of our results: k=500 and k=2000 for each of the two temperature field scenarios. Here, we consider the simulation to have reached a steady-state solution once the surface elevation changes by less than 0.02 m per year (which, with a timestep of 1.25 a, takes around several hundred iterations to achieve). Also, since the different model simulations result in firn slabs of different thicknesses, when plotting, we offset the model frames to ensure a common surface elevation, allowing for an easier comparison with the observed elevation profile. The steady-state surface elevations and BCO depths modeled are plotted in Fig 10b, c, respectively, along with the observed profiles. Overall, we find that the model (red and yellow lines) predicts higher rates of densification in the shear margins, where the shallowest part of the modeled firn columns are found to align with the maximum in horizontal shear strain rates. As expected, increasing kor shear-margin temperatures results in a softer firn that shortens the shear-margin firn column (raises the BCO depth). Comparing the surface elevation profiles (Fig. 10b), the model can at best explain half the depth of the shear-margin troughs, with the hot shear-margin scenario predicting the deepest troughs. Once a common vertical offset is set, the predicted surface https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press 16 Jon Arrizabalaga-Iriarte et al. stress and effective strain rate, . 𝜖E, should follow the power law . 𝜖E=A𝜎n E, corresponding to a Norton–Bailey creep potential. The flow rule then implies (Duva and Crow, 1994) . 𝝐= A 2𝜎n−1 E𝜕𝜎2 E 𝜕𝝈,(B2) where Ais the flow-rate factor. The derivatives needed for evaluating 𝜕𝜎2 E/𝜕𝝈 are 𝜕(𝝈 : 𝝈)/𝜕𝝈 = 2𝝈and 𝜕tr(𝝈)2/𝜕𝝈 = 2tr(𝝈)I, and the forward rheology is therefore . 𝝐=A𝜎n−1 E[a(𝝈− tr(𝝈) 3I)+2 3btr(𝝈) 3 I 3],(B3) 𝜎2 E=a 2(𝝈:𝝈− tr(𝝈)2 3)+b 3(tr(𝝈) 3)2,(B4) where a factor of 3(n+1)/2/2 has been absorbed into A. Notice that Glen’s incompressible flow law is recovered in the limit a=1 and b=0. The inverse rheology can be derived by vectorizing (B3) according to 𝒱(X)=[X11,X21,X31,X12,X22,X32,X13,X23,X33]T,(B5) giving 𝒱(. 𝝐)=A𝜎n−1 EP⋅𝒱(𝝈), (B6) where P=aI9+(2b 33−a 3)𝒱(I)⊗𝒱(I). (B7) Here, ⊗is the generalized outer product (Kronecker product), and tr(𝝈) = I: 𝝈 = 𝒱(I)⋅𝒱(𝝈)was used. The inverse rheology is then given by 𝝈 = 𝒱−1(P−1𝒱(. 𝝐)), or in tensorial form as written in Equation (4). https://doi.org/10.1017/jog.2025.6 Published online by Cambridge University Press