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
