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

Jou nal o Glaciology A icle Ci e his a icle: A izabalaga-I ia e J, Lejonagoi ia-Ga mendia L, H idbe g CS, G ins ed A, Ra hmann NM (2025) Fi n densi ica ion in wo dimensions: modeling he collapse o snow ca es and enhanced densi ica ion in ice-s eam shea ma gins. Jou nal o Glaciology 71, e61, 1–16. h ps:// doi.o g/10.1017/jog.2025.6 Recei ed: 4 Ap il 2024 Re ised: 6 Janua y 2025 Accep ed: 13 Janua y 2025 Keywo ds: A c ic glaciology; Ice physics; Ice heology; Ice s eams; Pola i n Co esponding au ho : Jon A izabalaga-I ia e; Email: jon.a izabalaga@bc3 esea ch.o g © The Au ho (s), 2025. Published by Camb idge Uni e si y P ess on behal o In e na ional Glaciological Socie y. This is an Open Access a icle, dis ibu ed unde he e ms o he C ea i e Commons A ibu ion licence (h p://c ea i ecommons.o g/licenses/ by/4.0), which pe mi s un es ic ed e-use, dis ibu ion and ep oduc ion, p o ided he o iginal a icle is p ope ly ci ed. camb idge.o g/jog Fi n densi ica ion in wo dimensions: modeling he collapse o snow ca es and enhanced densi ica ion in ice-s eam shea ma gins Jon A izabalaga-I ia e1,2,3, Lide Lejonagoi ia-Ga mendia2, Ch is ine S. H idbe g2, Aslak G ins ed2and Nicholas M. Ra hmann2 1Basque Cen e o Clima e Change (BC3), Leioa, Spain; 2Niels Boh Ins i u e, Uni e si y o Copenhagen, Copenhagen, Denma k and 3Facul y o Science and Technology, Uni e si y o he Basque Coun y (UPV/EHU), Leioa, Spain Abs ac Accu a e modeling o i n densi ica ion is necessa y o ice co e in e p e a ion and assessing he mass balance o glacie s and ice shee s. In his pape , we e isi he nonlinea - iscous i n heology in oduced by Gaglia dini and Meyssonnie (1997) ha allows mul idimensional i n densi ica ion p oblems o be posed, subjec o a bi a y s ess and empe a u e ields. Fi s , we ex end he cali- b a ion o he coe icien unc ions ha con ol i n comp essibili y and iscosi y o i e addi ional G eenlandic si es, showing ha he o iginal calib a ion is no uni e sally alid. Nex , we demon- s a e ha he ansien collapse o a G eenlandic i n unnel can be ep oduced in a c oss-sec ion model, bu ha anomalous wa m summe empe a u es du ing 2012–14 educe con idence in a emp s o independen ly alida e he heology. Finally, we show ha he heology can explain he inc eased densi ica ion a e and a ying bubble close-o dep h obse ed ac oss he shea ma gins o he No heas G eenland Ice S eam. Al hough we sugges mo e wo k is needed o cons ain he nea -su ace comp essibili y and iscosi y unc ions o he heology, ou esul s s eng hen he empi ical g ounding o he heology o u u e use, such as modeling ho izon al i n densi y a ia ions o e ice shee s o mass-loss es ima es o es ima ing ice-gas age di e ences in ice co es subjec o complex s ain his o ies. 1. In oduc ion Snow p ecipi a ed o e glacie s and ice shee s ans o ms in o ice in he uppe mos laye s o he ice column h ough i n densi ica ion. Du ing his p ocess, old snow is bu ied and comp essed as a esul o he o e bu den p essu e un il i eaches he densi y o pu e ice. Unde s anding and modeling he de ails o his p ocess is cen al o he in e p e a ion o ice co e eco ds and assessing he mass balance o glacie s and ice shee s based on sa elli e al ime y. In he case o ice co e eco ds, a good unde s anding o he ai - apping mechanism is key o in e p e ing gas eco ds. The ai p esen in i n po es becomes isola ed om he a mosphe e no un il deep in he i n column, a he so-called bubble close-o (BCO) dep h. This can c ea e 100- o 1000-yea di e ences in age be ween apped gases and he su ounding ice (Δage), which mus be conside ed o p ope ly synch onize eco ds (Schwande and o he s, 1997). In he case o mass-balance moni o ing om sa elli e al ime y, he p oblem is o ansla e changes in al i ude in o changes in mass. Sea-le el ise es ima es can be biased i he i n densi y ield (i.e. i n compac ion) is no co ec ly accoun ed o (Sø ensen and o he s, 2011, Lipo sky, 2022), and ecen e o s o include ho izon al densi y a ia ions o e he Wes An a c ic ice shee in la ge-scale modeling inds a co ec ion in olume abo e lo a ion o e 40 yea s o 10% (Schelpe and Gudmundsson, 2023). Fi n densi ica ion models a e an impo an ool ha can help es ima e he e ec o densi- ica ion on ice-co e clima ic eco ds and mass-loss unce ain ies. Ideally, densi ica ion models would be de i ed om i s p inciples, bu densi ica ion is a complex p ocess in which c ys- als ea ange and change hei shape and size. As a consequence, he ela i e impo ance o he di e en mechanisms d i ing densi ica ion is di ided in o se e al s ages depending on local condi ions ( o de ails see Cu ey and Pa e son, 2010, p. 16–22). To simpli y ma e s, la ge-scale models o i n columns he e o e end o be cons uc ed based on phenomenological a gumen s. He on and Langway (1980) p oposed a one-dimensional empi ically uned model o he i s and second s ages o densi ica ion o p edic ing dep h–densi y, dep h–load and dep h–age p o iles based exclusi ely on local empe a u e and accu- mula ion a e. Thei celeb a ed model has since become he benchma k o compa ing mo e sophis ica ed models wi h (e.g. Ba nola and o he s, 1991, Li and Zwally, h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess 2 Jon A izabalaga-I ia e e al. 2011). Howe e , ep oducing obse ed densi y p o iles a si es wi h di e se clima ic condi ions emains challenging, and some models gi e inconsis en esul s o he same bounda y condi ions (Lundin and o he s, 2017), sugges ing ha mo e wo k is wa an ed. Densi ica ion models ha a emp o expand on he physi- cal igo o p e ious wo k ha e also been p oposed, al hough imp o ed accu acy o e empi ically based models does no neces- sa ily ollow (Thompson-Munson and o he s, 2023). Fo example, Alley (1987) p oposed a simple model o po ous i n densi i- ca ion by g ain bounda y sliding, se e al models make use o a one-dimensional compac i e iscosi y (A naud and o he s, 2000, Mo is and Wingham, 2014, S e ens and o he s, 2023), and h ee- dimensional cons i u i e ela ions exis ha a gue o ce ain is- cosi y unc ions o i n (Salama in and o he s, 2009, Fou eau and o he s, 2024). Finally, we no e ha A he n and o he s (2010) p o- posed a Naba o–He ing c eep equa ion o po ous ma e ials ha is he basis o se e al o he models used o al ime y co ec- ions (Smi h and o he s, 2020), and SNOWPACK—a physically based land-su ace snow model o iginally de eloped o suppo a alanche wa ning (Ba el and Lehning, 2002, Lehning and o h- e s, 2002a,2002b)—has also been success ully adap ed o model nea -su ace i n densi ica ion (Keenan and o he s, 2021). 1.1. Beyond one-dimensional modeling A cen al p oblem is o ake in o accoun he la ge -scale ice low, as he e is e idence o accele a ed i n densi ica ion when ho izon al s ain a es a e nonze o (Ki chne and o he s, 1979, Alley and Ben ley, 1988, Ri e man and o he s, 2019). This phe- nomenon, called s ain so ening, was i s in oduced by Alley and Ben ley (1988); hey a gued ha du ing he second s age, densi- ica ion is p ima ily d i en by disloca ion c eep, which inc eases wi h he squa e o he e ec i e s ess. Hence, in egions whe e non- negligible ho izon al s ain a es a e supe imposed on he e ical i n compac ion, he e ec i e iscosi y is educed, which should, in u n, enhance densi ica ion. Mo i a ed by his, O aschewski and G ins ed (2022) ecen ly ex ended he He on and Langway (1980) model (HL) by in oduc- ing a scale ac o ha allows he inclusion o ho izon al s ain a es as a so o o cing pa ame e , oge he wi h he commonly used clima e a iables ( empe a u e and accumula ion a e). Speci ically, he densi ica ion a e de i ed om he classical HL model is mul i- plied by a scale ac o ha is a unc ion o he o he wise un esol ed s ain a e componen s. Al hough ex ensions based on He on and Langway (1980) a e gene ally compu a ionally inexpensi e, hey emain one-dimensional. Mul idimensional app oxima ions can be cons uc ed by ho izon ally joining i n-column models (e.g. O aschewski and G ins ed, 2022), bu o mula ing i n densi ica- ion in e ms o a h ee-dimensional cons i u i e s ess–s ain a e ela ionship is desi able o se e al easons: (i) he e ec o nonze o ho izon al s ain a es can be ep esen ed na u ally (wi hou he need o any ad hoc inclusion), (ii) wo- o h ee-dimensional i n compac ion p oblems can be mo e easily sol ed o complica ed bounda y condi ions and geome ies using, e.g. he ini e elemen me hod and (iii) Glen’s low law o solid ice may be ex ended o also cha ac e ize po ous i n, he eby allowing o a seamless nume ical ea men o he en i e i n–ice column. The semi-empi ical heology p oposed by Gaglia dini and Meyssonnie (1997)—hence o h GM97—is exac ly such a ield o mula ion. I is based on a gene al comp essible powe -law he- ology o po ous ma e ials, adap ed o i in si u measu ed densi y p o iles om Si e-2, G eenland (see Fig. 2), and cold oom de o - ma ion es s. The GM97 heology has p e iously been shown capa- ble o ep oducing age–dep h p o iles o moun ain glacie s (Lü hi and Funk, 2000, Zwinge and o he s, 2007, Gilbe and o he s, 2014, Licciulli and o he s, 2020) and p edic ing he e olu ion o snow ca es bu ied o e ime in An a c ica (B ondex and o he s, 2020). The GM97 heology depends on wo coe icien unc ions o densi y ha de e mine he local comp essibili y and iscosi y. The unc ional o m o hese, and associa ed calib a ion cons an s, ha e ecei ed li le a en ion in he li e a u e apa om he ini ial ea - men by Gaglia dini and Meyssonnie (1997) and subsequen ly by Zwinge and o he s (2007). The p oposed coe icien unc ions mainly disag ee o nea -su ace i n densi ies, which can ha e a la ge impac on modeled su ace eloci ies and densi ica ion a es (Gilbe and o he s, 2014). This has led o some disag eemen in he li e a u e abou which coe icien unc ions o use. Mo e wo k is he e o e wa an ed o make clea whe he he calib a ions p o- posed by Gaglia dini and Meyssonnie (1997) and Zwinge and o he s (2007) (i) hold o a wide ange o si es, whe e i n co es ha e since been d illed and (ii) hold in wo-dimensional se ings, whe e non i ial densi ica ion ields ha e been obse ed. I so, his would pu he GM97 heology on s onge empi ical g ounds and inc ease con idence in i s u u e use, which is he p ima y aim o ou wo k. In his pape , we es he lexibili y and use ulness o he GM97 heology o modeling i n compac ion beyond one-dimensional p oblems. Fi s , we e isi he o iginal calib a ion o he heology and assess he uni e sali y o i by ying o ep oduce i e addi- ional G eenlandic ice co e si es, whe e one-dimensional modeling is app op ia e. Then, we es he model’s accu acy and depen- dency on empe a u e by a emp ing o ep oduce he collapse o a nea -su ace snow ca e, cons uc ed a he NEEM ice co e si e in G eenland. Finally, we es whe he he heology can ep oduce he enhanced densi ica ion obse ed a he shea ma gins o he No heas G eenland Ice S eam (NEGIS), belie ed o be caused by s ain-so ening e ec s. In he ollowing, we begin by in oduc- ing he GM97 heology be o e conside ing he special one- and wo-dimensional cases he eo , needed o ou expe imen s. 2. Me hod Fo a comp essible ma e ial, mass conse a ion implies ha 𝜕𝜌 𝜕 +∇⋅(𝜌u)=0,(1) whe e uand 𝜌a e he eloci y and densi y ields, espec i ely. The eloci y ield o a S okes luid is go e ned by he momen um balance ∇⋅𝝈+𝜌g=0,(2) whe e 𝝈(. 𝝐)is he s ess enso , . 𝝐=(∇u+(∇u)T)/2 is he s ain a e enso and gis he g a i a ional accele a ion ec o . The GM97 heology ea s i n densi ica ion as a seconda y c eep p oblem in ol ing u,𝜌and he empe a u e (o en halpy) T, assuming negligible e ec s om p ima y c eep, snow me a- mo phism and he b i le ac u ing o snow. The heology ep e- sen s bo h i n and ice using a single, comp essible No on–Bailey powe -law luid ha con o ms o Glen’s incomp essible low law in he limi 𝜌→𝜌ice, whe e 𝜌ice =917 kg m−3is he densi y o glacie ice. The heology is o be unde s ood as an ex ension o Glen’s low law ha depends on he i s in a ian o he s ain a e enso , h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess Jou nal o Glaciology 3 (. 𝝐), in addi ion o he usual second in a ian , . 𝝐:. 𝝐. The heology was i s p oposed by Du a and C ow (1994) o gene al po ous ma e ials bu has since been adop ed in glaciology, and ecen ly wi h enewed in e es (B ondex and o he s, 2020, Licciulli and o h- e s, 2020). W i en in i s in e se o m (see Appendix B) equi ed by he momen um balance (2), he heology is 𝝈=A−1/n. 𝜖(1−n)/n E[1 a(. 𝝐− (. 𝝐) 3I)+3 2b (. 𝝐)I],(3) . 𝜖2 E=1 2a(. 𝝐:. 𝝐− (. 𝝐)2 3)+3 4b (. 𝝐)2,(4) whe e n=3 ollowing he abo e-men ioned adop ions o he he- ology, . 𝜖Eis he e ec i e s ain a e and aand ba e he coe icien unc ions con olling i n iscosi y and comp essibili y (elabo a ed on below). No e he e ha Glen’s law is eco e ed when a=1, b=0 and (. 𝝐)=0, which mus be ul illed in he limi 𝜌→𝜌ice. Since he low- a e ac o , A, depends on empe a u e, he e o- lu ion o he empe a u e ield (T) needs o be conside ed, oo: 𝜌c(𝜕T 𝜕 +u⋅∇T)=∇⋅(kT∇T)+𝝈:. 𝝐, (5) whe e c=c(T)and kT=kT(T,𝜌)a e he speci ic hea capac- i y and he mal conduc i i y o ice (see Appendix A o empi ical unc ional o ms used). In GM97, he low- a e ac o is modeled he usual way as an A henius ac i a ion unc ion o empe a u e (G e e and Bla e , 2009, G e e and o he s, 2014): A(T)=A0exp[−Q/(RT)], (6) whe e, ollowing Zwinge and o he s (2007), he canonical depen- dence on he empe a u e o he p e ac o A0and he ac i a ion ene gy o c eep, Q, a e p esumed: A0={3.985 ×10−13 s−1Pa−3 o T≤−10∘C 1.916 ×103s−1Pa−3 o T>−10∘C,(7) Q={60 kJ mol−1K−1 o T≤−10∘C 139 kJ mol−1K−1 o T>−10∘C.(8) The luidi y also depends on o he ac o s such as p essu e (Wee man, 1983), wa e con en (Ba nes and o he s, 1971, Du al, 1977), impu i ies (H o hold and o he s, 2012) and g ain sizes and o ien a ions (Du al, 1973, Shoji and Langway, 1988). The p essu e dependence is es ima ed o be ela i ely small, e en a hyd os a ic p essu es nea ice shee beds, and he e o e neglec ed (Rigsby, 1958). Rega ding wa e con en , empe a u es a e gene ally con- side ed o be low enough o neglec he p esence o liquid wa e . Impu i ies, such as dus , a e known o inc ease ice so ness by up o a ac o o wo, bu in he G eenland i n column conside ed he e— which consis s only o Holocene deposi ion— he dus loading is minimal o can a leas be assumed o be uni o m o i s o de . Finally, changes in he size and o ien a ion o g ains a e hough o be impo an in dynamic egions such as ice s eams and hei shea ma gins (e.g. Ge be and o he s, 2021), bu expanding (3)–(4) o also allow o iscous aniso opy is a daun ing ask and ou o scope o his wo k, so he e ec o de eloped c ys al ab ics in he i n cons i u e a heological unce ain y h oughou his wo k. Figu e 1. Coe icien unc ions a( ed) and b(blue) p oposed by Gaglia dini and Meyssonnie (1997) and Zwinge and o he s (2007) o n=3. 2.1. Coe icien unc ions Following Du a and C ow (1994), he coe icien unc ions aand ba e aken o depend exclusi ely on he ela i e densi y 𝜌= 𝜌 𝜌ice .(9) In his way, i aand ba e a ec ed by empe a u e, g ain sizes, e c., i is implici ly assumed ha such dependencies can be ac o ed ou and abso bed in o A. Indeed, his sepa a ion o dependencies is desi able since (3)–(4) can hen be easily made o educe o Glen’s law in he limi 𝜌→1, as no ed abo e. The coe icien unc ions aand bdepend in pa on a model o how s esses concen a e in he p esence o enclosed oids (ai ) (Du a and C ow, 1994) o densi ies abo e he c i ical alue 𝜌≥ 𝜌c i =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 pa on an empi ical scaling ela ion o densi ies 𝜌<𝜌c i ha has been ound o i cold oom expe imen s and densi ica ion measu ed a Si e 2, G eenland. Speci ically, wo such se s o aand b scalings ha e been p oposed: Gaglia dini and Meyssonnie (1997) sugges ed (Fig. 1, dashed lines) a1=(a0/b0)b1(12) b1=⎧ { { ⎨ { { ⎩ exp(451.63 𝜌2−474.34 𝜌+128.12) o 𝜌<0.5 exp(−17.15 𝜌+12.42) o 0.5≤ 𝜌≤ 𝜌c i b0 o 𝜌c i <𝜌 (13) h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess 4 Jon A izabalaga-I ia e e al. whe eas Zwinge and o he s (2007) sugges ed (Fig. 1, da k solid lines) a1={kaexp(−𝛾a( 𝜌− 𝜌s c)) o 𝜌≤ 𝜌c i a0 o 𝜌c i <𝜌 (14) b1={kbexp(−𝛾b( 𝜌− 𝜌s c)) o 𝜌≤ 𝜌c i b0 o 𝜌c i <𝜌 .(15) He e 𝜌s c =0.4 is he assumed ela i e densi y o su ace snow, con inui y a 𝜌= 𝜌c i implies he scaling exponen s a e 𝛾a=ln(ka/a0( 𝜌c i )) 𝜌c i − 𝜌s c and 𝛾b=ln(kb/b0( 𝜌c i )) 𝜌c i − 𝜌s c ,(16) and kaand kba e pa ame e s o be calib a ed agains obse a ions o expe imen s, sugges ed by Zwinge and o he s (2007) o be k≡ka=kb≃1000.(17) No e ha kaand kba e he alues o aand ba he su ace, 𝜌= 𝜌s c. As poin ed ou by B ondex and o he s (2020), kaand kbshould ideally be e-calib a ed on a case-by-case basis, al hough adop ing he cold oom and Si e 2 calib a ion p oposed by Zwinge and o h- e s (2007) (k≃1000) has been ound su icien in se e al cases (Gilbe and o he s, 2014, B ondex and o he s, 2020, Licciulli and o he s, 2020) and is ega ded as he canonical alue. In his wo k, we e isi he bes - i alue o kby conside ing addi ional one- and wo-dimensional model expe imen s ha a e e alua ed agains obse a ions. 2.2. Nume ics In he nume ical expe imen s ha ollow, we sol ed he coupled densi y, momen um and he mal p oblem using FEniCS (Logg and o he s, 2012), elying on New on’s me hod o sol e nonlin- ea i ies. Fo easons explained below, he ice s eam scena io is no he mally coupled, bu he mechanical p oblem is sol ed using he same me hod. The Jacobian o he esidual o ms ( equi ed o New on i e a ions) we e calcula ed using he uni ied o m lan- guage (Alnæs and o he s, 2015), used by FEniCS o speci y weak o ms o pa ial di e en ial equa ions (PDEs), which suppo s au oma ic symbolic di e en ia ion. All weak o ms a e p esen ed in Appendix A. Fo ou wo-dimensional expe imen s, meshes we e con- s uc ed using gmsh (Geuzaine and Remacle, 2009) and upda ed be ween ime s eps o e ol e he in e io and ex e io ee-su ace bounda y. 3. Valida ion o heology We conside ed h ee independen nume ical expe imen s o es he p oposed bes - i alue o he calib a ion cons an , k≃1000 (Zwinge and o he s, 2007). Fi s , we swep o e k o de e mine which alue can bes ep oduce obse ed i n densi y p o iles om a wide ange o G eenlandic ice co e d ill si es (case 1). Second, we sough a alue o k ha could bes ep oduce he measu ed collapse o a ench cons uc ed a NEEM (S e ensen, 2014) (case 2). Thi d, we a emp ed o ep oduce he obse ed enhanced densi ica ion o e he shea ma gins o NEGIS by a ying k(case 3). 3.1. Case 1: Rep oducing G eenlandic ice co es We conside ed a o al o six G eenlandic i n co es om he deep d illing si es a DYE-3, GRIP, NGRIP, NEEM, Si e 2 and Figu e 2. Si es o G eenlandic ice co es used o alida e he GM97 heology and de e mine k. Sa elli e-de i ed eloci ies om he MEaSUREs p og am a e shown in colo ed con ou s (Howa , 2020). Si e A (Fig. 2). Since i n empe a u es a e app oxima ely dep h- cons an a each si e (B éan and o he s (2017); Table 1), he p oblem educes o ha o sol ing o 𝜌and he e ical eloci y uz i he densi y and eloci y ields a e app oxima ed as ho izon ally homogeneous (simila o adi ional one-dimensional modeling, al hough his assump ion migh be less well jus i ied a dynamic si es). Unlike he ansien wo-dimensional p oblems conside ed in he ollowing, we sol ed di ec ly o he s eady s a es o (1)–(2) subjec o he bounda y condi ions 𝜌(z=0)=𝜌s c,(18) uz(z=0)=−. a,(19) uz(z=−H)=−. a𝜌s c 𝜌ice ,(20) whe e His he heigh o he modeled i n column, 𝜌s c is he su ace snow densi y and . ais he snow accumula ion a e a he si e (Table 1). The bounda y condi ion (20) assumes ha , in s eady s a e, he mass lux in o he i n column equals ha exi ing he h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess Jou nal o Glaciology 5 Table 1. Dep h-a e aged empe a u e, accumula ion a e, and su ace densi y o each si e conside ed. Accumula ion a es and su ace densi ies ollow om B éan and o he s (2017). Fo EGRIP, he accumula ion a e and su ace densi y da a a e e ie ed om Ka lsson and o he s (2020) and Schalle and o he s (2016), espec i ely. Si e Dep h a e age T(∘C) . a (kg m−2a−1)𝜌s c (kg m−3) Si e-2 −25.0*(Langway, 1970) 360 350.1 Si e-A (C ê e) −29.5(Clausen and o he s, 1988) 282 321.7 DYE-3 −21.0(Dahl-Jensen and o he s, 1998) 500 357.0 GRIP −31.7(Johnsen, 2003) 210 367.0 NGRIP −31.5(Dahl-Jensen and o he s, 2003) 175 299.9 NEEM −28.8(O si and o he s, 2017) 200 307.2 EGRIP −28.0*(Zuh and o he s, 2021) 130 290 *A e age su ace empe a u e was used due o he lack o a empe a u e p o ile. column, whe e Hmus be aken su icien ly la ge o gua an ee ha he bo om o he model domain is pu e ice, 𝜌(z=−H)=𝜌ice. We de e mined he bes - i alue o kby sweeping a wide ange o alues and calcula ing he esul ing model–obse a ion mis i . Speci ically, we conside ed k∈[1;1000]and calcula ed he oo - mean-squa e-e o (RMSE) be ween he modeled and obse ed densi y p o iles: RMSE =√ √ √ ⎷1 N N ∑ i=1(𝜌(zi)−𝜌obs(zi))2,(21) whe e Nis he o al numbe o disc e e dep h le els a which he mis i is e alua ed. Compu ing he pe cen age e o o each le el is also an op ion (weighing he misma ch wi h espec o he e e - ence alue i is de ia ing om), bu he esul s ob ained in his way a e i ually iden ical (da a no shown). Fig. 3 shows he RMSE mis i s calcula ed o he lowe -densi y sec ion o he i n column ( 𝜌<0.8), whe e kis expec ed o ha e he la ges in luence (see he abo e sec ion on he coe icien unc- ions). The esul s sugges a mis i minimum exis s somewhe e be ween k=100 and k=500 bu wi h disag eemen s be ween he si es; ha is, a global bes - i kmigh no be s ic ly applicable ac oss all si es (dis ega ding e ec s om seasonal empe a u e a ia ions no ea ed he e o simplici y). Fo e e ence, Fig. 4 shows he model pe o mance in he case o Si e 2 o di e en alues o k. No e ha he sensi i i y o he magni ude o kdec eases o ela- i e densi ies abo e 0.8. The unde es ima ion o densi y a dep h is p esen in all he si es modeled. 3.2. Case 2: T ench collapse a NEEM In an a emp o u he cha ac e ize he model dependence on k, we sea ched o al e na i e alida ion expe imen s ha in ol e nea -su ace (low-densi y) i n, since he choice o ka ec s com- p essibili y he mos he e. We ound ha he obse ed collapse o a cylind ical subsu ace unnel, cons uc ed a he NEEM ice co e si e, is well-sui ed o his pu pose. The unnel was cons uc ed o de e mine whe he subsu ace enches, buil ollowing he Camp Cen u y snow-blowing cas ing echnique (Cla k, 1965), bu using la ge in la able balloons ins ead, could se e as d illing and science enches o deep ice co e p ojec s. To judge i s easibili y, he unnel collapse a e was acked o 3 yea s a e i s cons uc ion. The ini ial and inal geome ies a e shown in Fig. 5 o e e ence, al hough he pic u es a e om he no he n end o he ench, while we modeled he sou he n end whe e measu emen s we e made mos consis en ly. Figu e 3. Model–obse a ion RMSE o G eenlandic densi y p o iles o k∈[1,1000]. All he cu es keep inc easing mono onically om k=1000 and onwa ds (no shown). Figu e 4. Modeled Si e 2 densi y p o iles o a ious k(colo ed lines) compa ed o obse a ions (B éan and o he s (2017); g ay do s). The e ec o a gi en kis mos e iden nea he su ace: he highe he alue o k, he as e he nea -su ace densi- ica ion. Fo 𝜌>0.8, he model gene ally unde es ima es densi ies a a gi en dep h o all k. The unnel was cons uc ed by blowing snow ou o a la ge ench and ins alling an in la able balloon in he open ca i y. A e in la ing i , snow was blown back in o he ench again, bu ying he balloon in a mo e closely packed (dense ) i n wi h a ypi- cal densi y o 𝜌 ench =550 kg m−3(B ondex and o he s, 2020, S e ensen, 2022). Finally, a e a couple o days, he balloon was de la ed and emo ed, p ese ing he cas ed, ci cula s uc u e o he unnel. We simpli y he ench collapse p oblem by conside ing a e - ical c oss-sec ion model, hus implici ly assuming ansla ional in a iance along he unnel (s ic ly speaking, only ele an h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess 6 Jon A izabalaga-I ia e e al. Figu e 5. NEEM balloon ench geome y 2 mon hs a e cons uc ion on 7 Augus 2012 (a) and 3 yea s la e on 27 May 2015 (b). No pic u e o he ench was a ailable o he model a ge yea , 2014, so 2015 is shown ins ead. Mo eo e , pic u es show he no he n end, whe eas we used he ench geome y o he sou he n end ha was mos consis en ly measu ed du ing he expe imen (no pic u es a ailable). Pic u es we e kindly p o ided by J. P. S e ensen and ep in ed wi h pe mission acco ding o he NEEM ice-co e p ojec media wai e . Figu e 6. Zoom-in o he ini ial geome y and densi y ield o he NEEM ench model. High-densi y snow was back illed in o he balloon ench, c ea ing a ha dened shell su ounding he unnel compa ed o he backg ound densi y ield. o e y long unnels). The ini ial geome y and densi y ield (Fig. 6) was cons uc ed ollowing he echnical epo by S e ensen (2014) and subsequen ly allowed o e ol e he mo- mechanically h oughou he du a ion o he es , app oxima ely 2 yea s long. No e ha he oo includes he epo ed ex a 1 m o back illed snow, bu ha he p esence o o he s uc u es in he ench (connec ing unnels, cabins, e c.) is no conside ed he e, which migh change esul s sligh ly and a e a sou ce o unce ain y. The ini ial backg ound densi y ield is aken o be equal o he smoo hed densi y p o ile o he NEEM co e (B éan and o he s, 2017). The snow accumula ion a e is se o he annual a e age o e he ench, which is double he e e ence measu emen a NEEM due o wind-d i en excess accumula ion (S e ensen, 2014). The ini ial empe a u e ield is aken o be uni o m and equal o he i n column a e age measu ed a NEEM, ⟨TNEEM⟩(Table 1). Howe e , du ing he 5 days ha i ook o cons uc he unnel, bo h he ench and he snow o be back illed we e exposed o highe - han- usual empe a u es (be ween −3 and −6∘C; S e ensen (2014)). Thus, he ini ial empe a u e o he i n ha su ounded he unnel was possibly up o 25∘C wa me han he i n-a e age-backg ound empe a u e (mos likely less, hough). In o de o assess he po en- ial unde es ima ion o he ini ial collapse a e, caused by p esc ib- ing oo cold i n, we he e o e also conside a second scena io in which all he back illed snow s a s ou wi h a empe a u e o −5∘C. This is jus an app oxima ion and does no in end o accu- a ely ep esen he ini ial empe a u e ield bu a he allows us o gauge he impac ha a much ho e ench migh ha e on ou esul s. The model domain has a wid h o L=20 m and a ime- e ol ing heigh p o ile h(x, ) ela i e o a ixed bo om bounda y loca ed a a dep h o H=30 m below he ini ial su ace. Bo h Land Ha e la ge enough o a oid bounda y condi ions a ec - ing he solu ion (judged by ial and e o ). The bo om bounda y is ixed by a ee slip condi ion and kep a he si e’s dep h- a e aged i n empe a u e. A he su ace, we he mally o ce he model by se ing he i n su ace empe a u e equal o measu e- men s om he local G eenland Clima ic Ne wo k wea he s a ion, TGCne ( )(Vandec ux and o he s, 2023). To es ima e he impac o no aking he su ace empe a u e e olu ion in o accoun , we also an he model by o cing i wi h a cons an i n su ace em- pe a u e, se equal o he a e age measu ed su ace empe a u e o e he modeled pe iod, ⟨TGCne ( )⟩. A pe iodic bounda y con- di ion is imposed on he le and igh sides o close he he mal p oblem. In summa y, he bounda y condi ions a e uz(x,z=−H)=0,(22) T(x,z=−H)=Ticeco e(z=−H), (23) h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess Jou nal o Glaciology 7 Figu e 7. Modeled e olu ion o he NEEM unnel om 10 June 2012 o 27 May 2014. (a) E olu ion o he unnel heigh o di e en alues o k(line colo s) and he mal scena ios (line s yles). (b) Co esponding unnel c oss-sec ions a he end o he simula ion. The la ge and smalle black cu es in panel b ep esen he ini ial and inal measu ed unnel dimensions, espec i ely. All simula ions a e he modynamically coupled bu di e in he su ace empe a u e bounda y condi ion. Expe imen s deno ed by solid and dash-do ed lines use measu ed ime-e ol ing su ace empe a u es, whe eas dashed lines deno e expe imen s whe e he a e age su ace empe a u e was imposed ( esul ing in p ac ically iso he mal condi ions). The ho ench scena io includes an ini ially ho e - han-a e age back illed ench (−5∘C), which aims o be mo e ep esen a i e o he eal ini ial condi ions. The h ee g ey do s in panel ashow measu ed unnel dimensions (S e ensen, 2014). The backg ound ed line in panel ashows he smoo hed hou ly empe a u e eco d om he si e’s wea he s a ion ha is pa o GC-Ne (Vandec ux and o he s, 2023). The o iginal eco d con ains posi i e empe a u es in he i s summe (da a no shown due o smoo hing). T(x,z=h, )={TGCne ( )i a iable clima e, ⟨TGCne ( )⟩=−27.7∘Ci a e age clima e. (24) Unlike he s eady-s a e model abo e, he op su ace bounda y is allowed o e ol e eely, as is he in e io unnel bounda y. Bo h bounda ies we e upda ed using he A bi a y Lag angian-Eule ian me hod p o ided by FEniCS, allowing mesh bounda y e ices o be displaced acco ding o he local p oduc be ween eloci y and ime-s ep size. In his me hod, su ace accumula ion is aken in o accoun by adding a cons an e ical eloci y componen , . a, o he su ace bounda y. No e ha , in e ec , su ace accumula ion g adually leads o an inc ease in he a e age heigh p o ile since no mass exi s he model domain; he domain bo om will, he e o e, also g adually end owa ds pu e ice. Since we a e only in e es ed in he ela i e de o ma ion o he ench geome y—and no he change in cen e - o - ench posi ion— his is o no conce n and does no a ec ou esul s. Due o la ge densi y g adien s causing nume ical ins abili ies, we assume ha he su ace accumula ion has a densi y o 𝜌 ench a he han 𝜌s c. The accumula ion a e was he e o e scaled by a ac o o 𝜌s c/𝜌 ench =0.56 o ensu e ha he o al mass deposi ed (and hence o e bu den load) is co ec . Al hough his assump ion educes he densi ica ion a e o he esh snow laye s, his does no ha e any di ec e ec on he mechanical e olu ion o he unnel below. Howe e , eplacing he snow wi h a hinne laye o dense (and, hus, less ai y) i n may cause he unnel o be mo e sensi i e o su ace empe a u e a ia ions. We sol e he ansien p oblem o he ou alues k= {100,500,1000,2000}using an Eule ime-s epping schemewi h a ime s ep o Δ =0.75 days. The esul s a e summa ized in Fig. 7, whe e he unnel heigh e olu ion and he inal c oss-sec ion a e shown o each case. Once again, he smalle kis, he s i e he i n and, hus, he less he unnel de o ms. Fo alues lowe han k=2000, he p edic ed de o ma ion is oo small o all he modeled scena ios. Fo highe alues, he unnel de o ms oo much, and he ceiling s a s o ca e inwa d (no shown he e). The ho e he su ounding i n is, he as e he unnel is ound o sh ink, as an icipa ed (dash-do ed compa ed o solid lines). The e is a delayed esponse in he collapse a es modeled o a i- a ions in su ace empe a u e, which is a consequence o he ime equi ed o he he mal pe u ba ion o each he dep h o he un- nel. Addi ionally, we also no e ha he sensi i i y o he unnel o su ace empe a u e a ia ions dec eases wi h ime as i becomes mo e isola ed below he accumula ing snow. 3.3. Case 3: NEGIS shea ma gin oughs The NEGIS is a cohe en s uc u e o as ice low, ini ia ing less han 150 km om he ice di ide and ex ending mo e han 600 km o he coas (G ins ed and o he s, 2022). Ele a ion maps ob ained om A c icDEM (Po e and o he s, 2018) show ha , compa ed o a e age su ace heigh a ia ions, he e a e ela i ely p onounced dep essions in he shea ma gins o NEGIS (Fig. 8), and ecen wo k by H idbe g and o he s (2020) has e i ied he dep hs o he oughs (dep essions) in A c icDEM using GPS s akes. The o igin o he shea ma gin oughs was e ealed in a ecen seismic su ey by Ri e man and o he s (2019) (su ey ansec is shown in Fig. 9 as a whi e line, and he su ace heigh in Fig. 10b as a blue line), inding ha he e is enhanced densi ica ion a he shea ma gins (blue con ou s in Fig. 10c). O aschewski and G ins ed (2022) la e showed ha hese dep essions migh pa ly be explained by s ain so ening, as he ho izon al eloci y g adien s a e la ge he e. h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess 8 Jon A izabalaga-I ia e e al. Figu e 8. A c icDEM su ace opog aphy ups eam o NEGIS, hillshaded using he 3D isualiza ion so wa e Blende . The EGRIP d ill si e is loca ed a he black ball. No e ha No h has been o a ed 135∘clockwise o make he iew o he ice-s eam ma gin mos clea , hence low is owa d he bo om pa o he igu e. Figu e 9. Sa elli e-de i ed su ace eloci ies a ound EGRIP om he MEaSUREs p o- g am (Howa , 2020) and ansec (whi e line) o he 79 geophones used o he seismic su ey pe o med by Ri e man and o he s (2019). Rep oducing he shea -ma gin oughs in a model o Ri e man’s ansec he e o e makes o a good es case o he GM97 heology (which accoun s o s ain so ening) and po en ially allows o ano he independen way o es ima e he bes alue o k. In addi ion o he obse ed shea -ma gin oughs (blue line in Fig. 10b), we also include he obse ed BCO dep h (i.e. dep h a which 𝜌=830 kg m−3) as a calib a ion a ge (whi e line in Fig. 10c). No e ha he abili y o model BCO dep hs migh be less sensi i e o he alue o kand mo e sensi i e o he unc ional o ms o aand b(densi ies a e la ge a he BCO). We model he ansec using a se up ha combines he app oaches om bo h he one-dimensional i n column and NEEM ench models. The yz domain conside ed is L=47.5 km wide and H=200 m all, including a 5 km bu e on bo h la e al sides o p e en he bounda ies om a ec ing he in e io solu ion (ex ension no shown in he igu es). The bounda y condi ions a e ee-slip on he le and igh bounda y: uy(y=0,z)=0,(25) uy(y=L,z)=0,(26) whe eas a non-ze o mass lux is imposed on he bo om bounda y ha balances he su ace-in eg a ed mass lux: uz(y,z=−H)=−. a𝜌s c 𝜌ice .(27) He e, His aken la ge enough o ensu e ha he densi y a he bo - om bounda y is ha o pu e ice. The densi y ield is subjec only o he su ace bounda y condi ion 𝜌(y,z=0)=𝜌s c.(28) The su ace densi y is conside ed o be 𝜌s c =343 kg m−3 (a e - age o he op 2 m o i n, Schalle and o he s (2016)) while he accumula ion a e is aken o be uni o m ac oss he ansec o simplici y, . a=130 kg m−2 y −1 (Ka lsson and o he s, 2020), despi e i may be up o a 20%highe in he shea ma gins due o d i snow apped by he oughs (Ri e man and o he s, 2019). The su - ace heigh p o ile h(y) was upda ed ollowing he usual kinema ic equa ion o ee su ace e olu ion: 𝜕h 𝜕 =. a+u(s) z−u(s) y𝜕h 𝜕y,(29) whe e . ais he su ace accumula ion a e, and u(s) yand u(s) za e he su ace eloci y componen s in he yz model domain. When e- meshing a e upda ing he su ace bounda y, he densi y ield was linea ly in e pola ed on o he new mesh. In he absence o any ho izon al a ia ion in bounda y condi- ions and ini ial s a e, he densi y ield would emain ho izon ally homogeneous in ime ( he column would e ol e solely due o g a i a ional compac ion). Howe e , addi ional ou -o -plane s ain a e componen s play an impo an ole in he NEGIS ansec : he shea ma gins expe ience la ge ho izon al x–yshea ( ed line in Fig. 10a) ha can cause signi ican s ain so ening, and he unk migh expe ience a sligh ex ending low (i any) in he along- low x-di ec ion (blue line in Fig. 10a). We he e o e supe - impose he obse ed . 𝜖xy(y)p o ile on ou y–zmodel domain by assuming i o be dep h in a ian , while neglec ing . 𝜖xx and . 𝜖yy as hey a e compa a i ely small. All s ain a es we e calcula ed based on sa elli e-de i ed eloci ies om he MEaSUREs p og am (Howa , 2020) and sligh ly smoo hed o a oid nume ical ins abili- ies. Finally, lacking in o ma ion abou he . 𝜖xz shea componen , we se i o ze o ollowing he shallow shel app oxima ion, commonly used o model ice s eam low. The e ec o he ou -o -model-plane componen . 𝜖xy on s ain so ening is included by ex ending he s ain a e in a ian s ha h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess Jou nal o Glaciology 9 Figu e 10. (a) Absolu e alue o he smoo hed ho izon al s ain- a e componen s along he NEGIS seismic ansec (solid lines), and he e ec i e ho izon al s ain a e . 𝜖e by O aschewski and G ins ed (2022), which includes he ups eam his o y, oo. (b) Modeled ( ed and yellow lines) and GPS-measu ed (blue line) su ace ele a ion anomaly p o ile along he seismic ansec (Ri e man and o he s, 2019). Solid and dashed lines ep esen iso he mal and ho shea -ma gin expe imen s, espec i ely. (c) Modeled BCO dep h p o iles (lines) plo ed on op o he obse ed densi y ield by Ri e man and o he s (2019). The whi e line shows he obse ed 𝜌=830 kg m−3BCO dep h con ou , and he iole line shows he BCO dep h modeled by O aschewski and G ins ed (OG22). Ve ical solid lines show he modeled shea -ma gin cen e poin s (deepes oughs), and dashed lines show he ho izon al ex en o he imposed 6∘C empe a u e anomaly in he shea ma gins. The along low dimension, x, is poin ing ou o he plane. en e he e ec i e iscosi y, . 𝜖E, acco ding o ( o be consis en wi h) hei h ee-dimensional de ini ions: (. 𝝐)= (. 𝝐2D), (30) . 𝝐:. 𝝐=. 𝝐2D :. 𝝐2D +2. 𝜖2 xy,(31) whe e . 𝝐2D is he wo-dimensional s ain a e enso in ou yz model. Fo he empe a u e ield ac oss he ansec , we conside wo di e en scena ios. In he i s scena io, we assume a uni o m em- pe a u e o T= −28∘C (and hus a uni o m low- a e ac o ) (Table 1). In he second scena io, we impose a 6∘C ele a ed em- pe a u e in he ice-s eam shea ma gins as epo ed by Holschuh and o he s (2019) (high-end es ima e, could be lowe ), pos u- la ed o be caused by s ain hea ing. Since ice is a good he mal insula o , he 6∘C anomaly is speci ied as an ab up bu dep h- cons an s ep unc ion c ossing in o he shea ma gins. Al hough adding a he mal coupling o his p oblem as well (i.e. e ol ing he empe a u e ield) would inc ease he physical ealism, he ou -o - plane hea low is poo ly known, making he he mal p oblem no well-posed. On he o he hand, by imposing he abo e-men ioned s eady empe a u e ields, he e ec o dynamic s ain so ening (due o nonlinea iscosi y) is made clea e , as i s e ec can be unde s ood in isola ion. Fo b e i y, we epo only on ou ansien simula ions ha a e chosen o be ep esen a i e o ou esul s: k=500 and k=2000 o each o he wo empe a u e ield scena ios. He e, we con- side he simula ion o ha e eached a s eady-s a e solu ion once he su ace ele a ion changes by less han 0.02 m pe yea (which, wi h a imes ep o 1.25 a, akes a ound se e al hund ed i e a- ions o achie e). Also, since he di e en model simula ions esul in i n slabs o di e en hicknesses, when plo ing, we o se he model ames o ensu e a common su ace ele a ion, allow- ing o an easie compa ison wi h he obse ed ele a ion p o- ile. The s eady-s a e su ace ele a ions and BCO dep hs modeled a e plo ed in Fig 10b, c, espec i ely, along wi h he obse ed p o iles. O e all, we ind ha he model ( ed and yellow lines) p e- dic s highe a es o densi ica ion in he shea ma gins, whe e he shallowes pa o he modeled i n columns a e ound o align wi h he maximum in ho izon al shea s ain a es. As expec ed, inc easing ko shea -ma gin empe a u es esul s in a so e i n ha sho ens he shea -ma gin i n column ( aises he BCO dep h). Compa ing he su ace ele a ion p o iles (Fig. 10b), he model can a bes explain hal he dep h o he shea -ma gin oughs, wi h he ho shea -ma gin scena io p edic ing he deepes oughs. Once a common e ical o se is se , he p edic ed su ace h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess 16 Jon A izabalaga-I ia e e al. s ess and e ec i e s ain a e, . 𝜖E, should ollow he powe law . 𝜖E=A𝜎n E, co e- sponding o a No on–Bailey c eep po en ial. The low ule hen implies (Du a and C ow, 1994) . 𝝐= A 2𝜎n−1 E𝜕𝜎2 E 𝜕𝝈,(B2) whe e Ais he low- a e ac o . The de i a i es needed o e alua ing 𝜕𝜎2 E/𝜕𝝈 a e 𝜕(𝝈 : 𝝈)/𝜕𝝈 = 2𝝈and 𝜕 (𝝈)2/𝜕𝝈 = 2 (𝝈)I, and he o wa d heology is he e o e . 𝝐=A𝜎n−1 E[a(𝝈− (𝝈) 3I)+2 3b (𝝈) 3 I 3],(B3) 𝜎2 E=a 2(𝝈:𝝈− (𝝈)2 3)+b 3( (𝝈) 3)2,(B4) whe e a ac o o 3(n+1)/2/2 has been abso bed in o A. No ice ha Glen’s incomp essible low law is eco e ed in he limi a=1 and b=0. The in e se heology can be de i ed by ec o izing (B3) acco ding o 𝒱(X)=[X11,X21,X31,X12,X22,X32,X13,X23,X33]T,(B5) gi ing 𝒱(. 𝝐)=A𝜎n−1 EP⋅𝒱(𝝈), (B6) whe e P=aI9+(2b 33−a 3)𝒱(I)⊗𝒱(I). (B7) He e, ⊗is he gene alized ou e p oduc (K onecke p oduc ), and (𝝈) = I: 𝝈 = 𝒱(I)⋅𝒱(𝝈)was used. The in e se heology is hen gi en by 𝝈 = 𝒱−1(P−1𝒱(. 𝝐)), o in enso ial o m as w i en in Equa ion (4). h ps://doi.o g/10.1017/jog.2025.6 Published online by Camb idge Uni e si y P ess