A cut finite element method for non-Newtonian free surface flows in 2D: application to glacier modelling

Contents lists available atScienceDirect

Journal

of

Computational

Physics:

X

www.elsevier.com/locate/jcpx

A cut

finite

element

method

for

non-Newtonian

free

surface

flows

in

2D

- application

to

glacier

modelling

Josefin Ahlkrona

a,b,

∗

,

Daniel Elfverson

c

a_{DepartmentofMathematics,StockholmUniversity,Stockholm,Sweden} b_{Swedishe-ScienceResearchCentre(SeRC),Stockholm,Sweden}

c_{Department ofMathematicsandMathematicalStatistics,UmeåUniversity,Umeå,Sweden}

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received27August2020

Receivedinrevisedform25February2021 Accepted4April2021

Availableonline15April2021

Keywords: Icesheetmodelling CutFEM

Freeboundaryproblems Non-Newtonianflow

Unfittedfiniteelementmethods Sharpinterfacemethods

Inicesheetandglaciermodelling,theFiniteElementMethodisrapidlygainingpopularity. However, constructing and updating meshes for ice sheets and glaciersis a non-trivial and computationally demanding task due to their thin, irregular, and time dependent geometry. In this paper we introduce a novel approach to ice dynamicscomputations based ontheunfitted FiniteElement MethodCutFEM,whichletsthe domainboundary cutthroughelements.ByemployingCutFEM,complexmeshingandremeshingisavoided as the glaciercan beimmersed inasimplebackground meshwithoutlossofaccuracy. The iceismodelledas anon-Newtonian,shear-thinningfluidobeying thep-Stokes(full Stokes)equationswiththeiceatmosphereinterfaceasamovingfreesurface.ANavierslip boundaryconditionappliesattheglacierbaseallowingbothbedrockandsubglaciallakes toberepresented.WithintheCutFEMframeworkwedevelopastrategyforhandling non-linearviscositiesand thin domainsandshow howglacierdeformationcanbe modelled usingalevelsetfunction.Innumericalexperimentsweshowthattheexpectedorderof accuracyisachievedandthatthemethodisrobustwithrespecttopenaltyparameters.As anapplication wecomputethevelocityfieldoftheSwissmountainglacierHautGlacier d’Arolla in2Dwith and withoutanunderlyingsubglacial lake,and simulatetheglacier deformationfrom year1930to 1932, withand withoutsurfaceaccumulationand basal melt.

©2021TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCC BYlicense(http://creativecommons.org/licenses/by/4.0/).

1. Introduction

Asthe atmosphereandoceaniswarming,icesheetsandglaciersmelt.Themeltwater causesglobalmeansealevelto rise, andhas thepotential toalter ocean circulation [1,2]. Icebehaves as a very viscous,non-Newtonian, shear-thinning fluid,andinordertopredictfuturemeltratesitiscrucialtobeabletoaccuratelysimulatehowitflowsanddeforms.Ice sheetandglaciermodellingisthusaveryimportantmovingboundary,non-linear,fluidmechanicsproblem.

Due to the complexandchanging geometries involved, theFinite Element Method(FEM) isgaining popularityin ice modellingsince itallows forunstructuredmeshes[3–6].Icesheet andglaciermodels employ traditionalfitted FEM,that is,thecomputationalmeshisfittedtothemodeldomainwithnodesplacedontheboundary.However,thereisalimitto justhowcomplexmovingdomainsfittedFEMcanhandle.Requiringmeshnodestoalignwiththeboundaryinvolves

non-*

Correspondingauthor.

E-mailaddress:ahlkrona@math.su.se(J. Ahlkrona).

https://doi.org/10.1016/j.jcpx.2021.100090

2590-0552/©2021TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

trivialmeshconstruction,andasthedomainmoveseithercostlyremeshingorproblemswithlow qualitymeshelements ifthemesh issimplydistorted.In icemodellingthe initialmesh isusually constructedby creatinga2Dtriangulationin the horizontalplane andextrudingit inthe verticaldirectioninto anumber oflayers. Becauseglaciers are thinnerthan theyarewide,theresultingprismaticelementsbecomeveryflat.Atlandbasedicemarginsthethicknessofelementseven goestowardszeroorsomeartificialminimumthickness,andwhentheice/atmosphereinterfacemoves,themeshisusually distortedto fitthe boundary,which cancomprise elementquality evenmore. Atthe coast,the extrudedmeshapproach rendersicecliffsatthecoastexactlyvertical,whichisnotalwaystrueinreality.Furthermore,whenthemarginsoftheice sheet movethemeshmayhavetobecompletelyreconstructed,whichisacostlyprocedure.Icesheet marginsmaymove suddenly, e.gwhenicebreaksoffintoicebergs,andinlongsimulations therecanevenbe topologicalchanges ifthe ice meltsintoseparateislands.

Inrecentyears,morerobustalternativestotraditionalfittedFEMhasbeendeveloped;unfittedsharpinterfacemethods [7,8]. Unfitted methods does not rely on placing mesh nodes on the boundary, but rather immerses the domain onto a background mesh,letting the boundary cut through elements. Apartfromavoiding generation andrepeated update of complexmeshes, unfittedmethods alsoallowsfortopologicalchanges.Inthispaperwe extendtheunfittedFEMmethod calledCutFEM[7],sothatitcanbeusedforicemodelling.

CutFEM relieson implementingboundary conditionsweakly using Nitsche’s method[9].Using only Nitsche’s method withan unfitted meshwouldtypically resultinlosinghalf an orderofconvergence[10]. CutFEMinstead keepsaccuracy byemploying quadratureonlyoverthephysicaldomain,combinedwithextrastabilisationterms.CutFEMhassuccessfully beenappliedtoStokesandNavier-Stokesflow[11,12],butmodellingicerequiresustoextendthemethodtothestrongly non-linearp-Stokes(orfullStokes)equationswhichdescribesnon-Newtonianflow,andtoadaptthemethodtothethinness andsharpcornersoficedomains.Inthispaperwe describehowthiscanbedoneandprovetheaccuracyofourmethod throughnumericalexperimentsinMATLAB.Weimplementalevel-setmethodformovingtheice/atmosphereinterfaceand show howtorepresentvariableslidingatthebaseoftheice.Nitsche’s methodandthelevel-setmethodhassuccessfully beenappliedtoicedynamicsmodellingine.g.[13,14] forpartsoftheboundary(thecalvingfrontandgroundingline),but outsidetheframeworkofCutFEM.TheadvantageofusingCutFEMisthatoptimalaccuracyiskept.

The paperisorganisedasfollows.InSection 2wedescribe thenon-linear equationsgoverningiceflowtogether with the equation that describes the moving boundary. In Section 3 theCutFEM method isoutlined, including how the non-linear materialproperties,thin domainsand basalsliding conditionsare handled.Section 4 explainsthe implementation of thelevel setmethodandhow glaciermargins arehandled astheicemoves. In Section 5we perform fournumerical experiments:oneunitsquaretesttoensurethattheorderofaccuracyiskeptandthattheaccuracyisrobustwithrespect topenaltyparameters,andsixexperimentsmodellingtheSwissmountainglacierHautGlacierd’Arolla,whichispartofthe icesheetmodellingbenchmarkexperimentISMIP-HOM[15].

2. Governingequationsandboundaryconditions

Glacierandicesheet flowisgovernedby anon-linearversionoftheStokes equationscalledthep-Stokes equationsor fullStokes equations,

−∇ ·

S

(

Du

)

+ ∇

p

=

ρ

g

∇ ·

u

=

0



in

.

(1)

Hereu

= (

ux

,

uy

)

isthevelocityfieldinatwodimensionalCartesiancoordinatesystem

(

x

,

y

)

, p isthepressure,



∈ R

2 is an open,two-dimensionalboundeddomain,and

ρ

g istheforce ofgravity. ThestresstensorisdenotedbyS anddepends onthestrainratetensorDu

= (∇

u

+ (∇

u

)

T

)/

2 as

S

(

Du

)

=

2

μ

Du

.

(2)

Sinceiceisanon-Newtonian,shear-thinningmaterial,theviscosity

μ

dependsontheFrobeniusnormofthestrainrate

μ

=

1 2A 1−p

₍

_

2 crit

+

1 2

||

Du

||

2 F

)

p−2 2

_.

₍₃₎

Thematerialparameter

p

liesintheinterval

(

1

,

2

)

forshear-thinningfluids,meaningtheviscosityisverylargewherethe strain rate islow. Asmall

p

corresponds to a strong non-linearitywhile

p

=

2 modelsa linear, Newtonian fluid. Forice

p

=

4

/

3 isthe standard choice, whichin glaciologyterminologycorresponds tothe Glenparameter beingnGlen

=

3.The smallregularisationparameter



critisincludedtoavoidnumericalproblemsincaseswhentheFrobeniusnormofthestrain ratetensor

||

Du

||

2_F approacheszero.Theratefactor A

=

A

(

T

)

modelshowiceviscositydependsonthetemperatureT .We considerisothermalconditionandsothat A isaconstant.

The natural velocity and pressure spaces for weak solutions ofthe full Stokes equations are

[

W1,p

_()

_]

2 _{and L}p

_()

respectively,where Lν denotesaLebesguespaceand W1,ν aSobolevspaceoforderone.Theparameter

p

 relatesto

p

as

p



= p/(p

−

1

)

,seee.g. [16,17].Notethatinalinearsetting

p

= p



=

2 andthespacesareconsistentwiththoseassociated withthelinearStokesequations.

Fig. 1. Across-sectionofthemountainglacierHautGlacierd’Arolla.Theoutwardspointingnormalisdenotedbyn andthetangentialvectorbyt.Thebase brestspartlyonbedrock(brown)andpartlyonsubglacialwater(blue).Thesurfacesevolvesintimeduetoicedeformationandprecipitation/melting.

The glacierboundary



= ∂

is divided into two parts,



s and



b, denoting the interface between the ice and the atmosphere(icesurface)andtheinterface betweenthegroundandtheicerespectively,seeFig.1.Theicesurface



s isa freeboundary,whichinglaciologyistypicallydescribedbyaheightfunction z

=

h

(

x

,

t

)

indicatingtheelevationof



s ata certainhorizontalpositionx atatimet.Inthispaperwe willinsteadusealevel-setformulation,whichunliketheheight functionformulationdoesnotrequirethatthereisauniquez foreachx,norexcludestopologicalchanges.Thepositionof thesurfaceisthusindicatedbyalevelsetfunction

(

x

,

z

,

t

)

whichiszeroatthesurface



s.Theevolutionofthesurfaceis givenbytheevolutionequation

∂



∂

t

+

u

· ∇ =

as

,

(4)

whereas isan accumulation/ablationfunction depending ontheatmosphericconditionssuch ase.g. snow.Ascustomary inglaciologytheaccumulation/ablationisaddedonlyintheverticaldirection(sinceprecipitationandablationcomesfrom above),i.e.as

= −

a⊥s

∇

,a⊥s beingmateriaaddedorremovedinthez-coordinatedirection.

TheboundaryconditionforthefullStokesequationsatthefreesurfaceis

(

−

pI

+

2

μ

Du

)

n

=

0 on



s

,

(5)

wheren istheoutwardspointingunit normalofthedomainboundaryandI istheidentitymatrix.Atthebase



btheice caneitherbefrozentothegroundoritcanbeslidingonsedimentsorwatersothataNaviersliptypeconditionapplies

Pn

(

u

−

gd

)

=

0 on



b

,

(6)

Pt

(



(

2

μ

Dun

−

d

)

+

μ

(

u

−

gd

))

=

0 on



b

.

(7) Here the projection matrices Pn

_,

_Pt

_{∈ R}

2×2 _{defines the projection to the normaland tangential plane and are given as}

Pn

=

n

⊗

n andPt

=

I

−

n

⊗

n,where

⊗

denotesan outer product. Weset d

=

0, meaningthere iszerotraction atthe base, andwillfromnowonomitd from theequations.Thetermgd

= (

0

,

−

a⊥_b

)

models meltorrefreezingatthebaseof theglacier.Thesliplength



determines theamountofslipinthetangentialdirection. For



=

0 equations(6)-(7) reduce to ahomogeneous Dirichletcondition, validforthecasewhen theiceisfrozen tounderlyingbase.When



→ ∞

afree slip condition is retrieved, which is valid when the ice is slipping over water. An intermediate slip length



results in a mixed Robincondition, which describessituations when the iceis not frozen atthe baseor is resting ondeformable debris/sediments.TheNavierslipboundaryconditionsofequation(6)-(7) wereformulatedinaCutFEMframeworkby[18] forthelinearStokesequations(

p

=

2).

3. CutFEMformulation

3.1. Computationalmeshandfiniteelementspaces

The physical domain



is immersed ina larger a rectangular domain



0,i.e.



⊂ 

0

⊂ R

2, see Fig. 2a. This larger domainiseasytopartitionintoastructuredsocalledbackgroundmesh

T

0,hconsistingofrectanglesK ,whereh denotesthe cell diameter.FortheArollaglacier, thebackgroundmeshes,

T

0,h are slightlyanisotropic toaccommodate thethinArolla glacier,unlikethebackgroundmeshesusedsofarintheframeworkofCutFEM.Elementsarethuseighttimeslongerinthe x-direction than inthe y-direction,reflecting thethinnessoftheglacier(see Fig.1). Thediscrete solutionwillbe sought

Fig. 2. AmeshfortheArollaglacier.Thecyanpartmarksthephysicaldomain.Thethinblacklinesmarksthemeshfacesforthebackgroundmesh constructedon0andthemagentaandcyantogethermarkstheactivepartofthemesh,Th().Notethatthepictureisstretchedintheverticaldirection, sothatinrealitytheelementsareslightlyanisotropic.

onlyontheactivepartofthemesh,definedasallelementswhichareintersectedbythedomain



,i.e.

T

h

()

= {

K

∈

T

0,h

:

K

∩ (

t

)

= ∅}.

(8)

Theactive meshismarkedcyanandmagentainFig.2a.Wecalltheunionoftheelementsintheactivepartofthemesh theactivedomain



h

(

t

)

= ∪K

∈ThK .Thefacetsoftheelementsin

T

h

()

entersintheformulationofstabilisation termsand wethereforeintroducethesomenotationsforthese.Allfacetsbelongingtoelementsin

T

h

()

aredenotedby

F

h

()

and all interior facets,i.efacets thatare sharedby two elements that areboth in

T

h

()

are denoted by

F

i,h

()

.The interior facetswillbeusedforstabilisingthepressure(i.e.forinf-supstabilisation).Theelementsandfacets

T

h

()

and

F

h

()

which are intersected,or“cut”,by theboundary



areofextraimportance,asitisonthesefacetsthatextrastabilisationterms areneededfortheCutFEMformulation.Theseelementsandfacetsaredefinedas

T

h

()

= {

K

∈

T

h

()

:

K

∩  = ∅},

(9)

F

h

()

= {

F

∈

F

i,h

()

:

K+F

∩  =

0 or K−F

∩  =

0

},

(10)

where K+_F andK−_F arethetwoelementsthatsharethefacet F .InFig.2bfacetsbelongingto

F

h

()

aremarkedyellow. Onthebackgroundmesh



0 weintroduceafiniteelementspaceofbi-linear,continuousfunctions

X0,h

:=



ω

∈

C

(

0

)

;

ω

K

◦

FK

∈ Q

1

( ˆ

K

),

K

∈

T

0,h



,

(11)

where the space

Q

1

( ˆ

K

)

:=

span

{ˆ

x

,

y

ˆ

,

ˆ

x

ˆ

y

,

(

x

ˆ

,

y

ˆ

)

∈ ˆ

K

}

is defined on the reference element K

ˆ

:= (

0

,

1

)

d. The mapping FK mapsfromthereferenceelementtoan element K ,andwillforthenumericalexperimentsontheArollaglaciergeometry representaslightcompressionofthegridintheverticaldirection.Thefiniteelementspaceontheactivedomainisdefined astherestrictionof X0,h onto



h

Xh

=

X0,h

|

h

.

(12)

Thisspaceisusedtoapproximateboththevelocityandpressuresolution

X

h

= [

Xh

]

2

,

Q

h

=

Xh

,

(13)

thatis,weemployequalorderbilinearelementsontheactivemesh.

3.2. Discreteproblem

Thediscretized problemnowreads:Findvelocityandpressure

(

u

,

ph

)

∈ X

h

×

Q

hsothat

ah

(

uh

,

vh

)

+

b

(

ph

,

vh

)

−

b

(

qh

,

uh

)

+

sp

(

ph

,

qh

)

(14)

whereah stems fromthestresstensor in(1), thetwo bh stemsfromthepressure gradientandfromthedivergencefree constraint, sh isan inf-supstabilisation neededaswe employequalorderelements, gu,and gp aretheextra stabilisation termscalledghostpenaltyterms neededfortheCutFEMformulation, andlh representstherighthandside of(1).We will nowdescribeeachofthesetermsindetail.

3.2.1. Nitscheterms

Following[18] theformah

(

uh

,

vh

)

containstheviscoustermtogetherwithseveralextraboundarytermsthatarerequired forimposingtheDirichletorNavier-slipboundaryconditions(6)-(7) weaklyonaboundarythatcutsthemesh.

ah

(

uh

,

vh

)

:= (

2

μ

Duh

,

Dvh

)



− 

2

μ

Duhn

,

vh



b

− 

P n_u h

,

2

μ

Dvhn



b

+ 

2

μ

γ

n_h  uh

·

n

,

vh

·

n



b

+ 





+

γ

t_h  Pt

(

2

μ

Duhn

),

vh



b

+ 

μ



+

γ

t_h  uh

·

t

,

vh

·

t



b

− 

γ

th



+

γ

t_h  Pt

(

2

μ

Duhn

),

2Dvhn



b

− 

μγ

th



+

γ

t_h  Ptuh

,

2Dvhn



b

.

(15) Thesecondtermcomesfromthepartialintegrationassociatedwithwriting (1) onweakform.Notethatunlikeforstandard fittedFEM,thistermisnon-zeroevenforDirichletconditionssincev isnon-zeroon



b.Thetermmakestheform unsym-metric,whichwithNitsche’smethodisremediedbyaddingextratermstosymmetrizetheproblem.Forno-slipconditions (



→

0) thethird andlast term together constitutesthesesymmetrisation terms.The fourthand sixth termare penalty terms that impose the Dirichlet condition, and

γ

t _and

_γ

n _{are Nitsche penalty parameters in the tangential and normal} direction.Thelowertheseparametersare,themoreemphasisisputonimposingtheboundaryconditions,versusfulfilling theequations.For



→ ∞

we retrievethe freeslipcondition.Inthiscaseitisthesecond tolasttermthatisthepenalty termimposingthe boundarycondition.Fora moredetaileddiscussionoftheseterms,werefer to[18]. Notethatforice, theviscosityineachofthesetermsisvariableintimeandspace,unlikeintheformulationin[18].Themeasureh differs from[18] wherethecellsizeh isusedinstead.Sinceforanisotropicelementsitmattersiftheboundary



cutsanelement alongthelongsideorshortside,wedefineh

=

hxhy

/(

hxt,x

+

hyt,y

)

,wheret istheaveragetangentialto



withinan element.FortheArollaglacier,theboundary mostlycutsthecells alongthelongedge ofelements,sothat h

≈

hy.This choicecanbemotivatedbyrevisingthetrace andinverseinequalitiesinvolvedin theCutFEMstability proofsandwillbe presentedthoroughlyinaseparatepublication.

Forthepressure/divergenceconstraintformbh oneextraboundarytermenters,stemmingfromthepartialintegration

bh

(

ph

,

vh

)

:=

b

(

ph

,

vh

)

+ 

p

,

v

·

n



b

.

(16)

Finally,giventhatd iszeroin(6)-(7),therighthandsideformis

lh

(

u

,

v

)

:= (

f

,

v

)



− 

gd

·

n

, (

2

μ

Dvhn

)

·

n



b

+ .

2

μ

γ

n_h  gd

·

n

,

vh

·

n



b (17)

− 

gd

·

n

,

qh



b

+ 

μ



+

γ

t_h  gd

·

t

,

vh

·

t



b

− 

μγ

t_h 



+

γ

t_h  Ptgd

,

2Dvhn



b

Notethateventhough thefiniteelementspacesaredefinedontheentireactive domain



h,theintegrationsin(15)-(17) areonlyover



foraccuracyreasons.ThecontributiontothefinalFEMmatricesofthecutelementsmaythusbeverysmall (seeFig.2b),whichwouldleadtobadconditionnumbersofthesematricesandstabilityissueswereitnotforgu,andgp.

3.2.2. Stabilisation terms

Theghostpenaltystabilisationtermsguandgpthatareaddedtohandleissuesrelatedtosmallcutelementsaredefined as gu

(

uh

,

vh

)

= β

u



F∈Fh() hn





nf

· ∇

uh



,



nF

· ∇

vh





F

,

(18) gp

(

ph

,

qh

)

= β

p



F∈Fh() h3 n

μ





nF

· ∇

ph

,



nF

· ∇

ph





F

,

(19)

where hn and ht is the length of the side normal and tangential to the facet F (see Fig. 3), nF is the normal vector to the facet,

β

u and

β

p are user definedparameters and

φ

:= (φ

+F

− φ

−F

)

denotes the jump of a quantify (

φ

), where

φ

±

(

x

,

y

)

:=

limt→0+

φ ((

x

,

y

)

±

tnF). The ghost penalties penalisejumps in pressure and velocity gradients across facets belonging to

F

h

()

. These stabilisationterms are an anisotropic version ofthe isotropic ghost penalty terms in [11], as weneedslightlyanisotropicbackgroundmeshestoaccommodatethethinArollaglacier.Theformulationcanbefoundby following thestandardCutFEM stabilityproof andtrackingthe lengthofthe normalandtangentialpartin thetrace and inverseinequalitiesseparately.

Fig. 3. The length of the sides of each cell is denoted by htif it is parallel to the facet F , and hnif it is perpendicular.

Thetermsp

(

uh

,

vh

)

in(14) isastabilisation termofContinuousInteriorPenalty(CIP)type[19]

sp

(

ph

,

qh

)

=

α



F∈Fh()

h2_xhn

μ





nF

· ∇

ph

,



nF

· ∇

ph





F

,

(20)

where

α

isauserdefinedparameter.Thisstabilisationisneededsincetheequalorderspaces

X

h and

Q

h arewell-known toviolate theinf-supcondition [20].Alsoherewehavemodified theformulationtosuitanisotropicmeshes. Theformof thestabilisation isconsistentwiththeanisotropicstabilisation of[21].Anindepthanalysisoftheanisotropicstabilisation terms(18) and(20) willbepresentedinafollowuppublication.

3.3. Handlingofthenon-linearity

Toresolvethe non-linearityofthe fullStokesequationswe useaPicard iteration,seeAlgorithm1.We choseaPicard iterationover aNewtoniteration,sincethefocusofthispaperisnot tomeasureefficiency.Newtoniterationsareusedin icemodellingbutrequiressomeextracareduetotheviscositysingularityinherenttothefullStokesequations.

Ineachiteration,theviscosity

μ

isupdatedintheviscoustermah,theghostpenaltyandtheinf-supstabilisationterm. Note that

μ

enters inall of theterms ofah.To ourknowledge CutFEMhas not beenused fora viscositythat varies in spacean time. We willshow by numericalexperimentsthat simplyupdating the viscosityinall termswhere itappears gives an accurate solution. The resulting linear system is solved using MATLAB’s inbuilt backslash operator

mldivide

and the computed velocity is then used in order to update the free surface in the manner described in the next sec-tion.

Algorithm1Generalsolutionprocedure.Foreachtimestepk andnon-lineariterationn alinearsystemissolved. 1: Setinitialconditionforvelocityu0

0,pressurep 0 0. 2: for k=1,...kmaxdo

3: n=1

4: while change>tolandn<50 do 5: Computeviscosityμk

n=μ(ukn)

6: Assemblelinearsystemfromproblem(14) usingμk n 7: Solvesystemforvelocityukn+1andpressurep

k n+1. 8: calculatechange= uk n+1−ukn2/ukn+12 9: n=n+1 10: end while 11: Useuk

ntoupdatethefreesurfaceks+1 12: UpdatetheactivemeshTh() 13: k=k+1

14: end for

4. Timeevolutionoftheicesurface

Inordertotrackthepositionofthefreesurface



sweuselevelsetfunction



: 

0

→ R

.Itisasigneddistancefunction, sothat

⎧

⎨

⎩

 <

0 if

(

x

,

y

)

∈ 

₊

,



=

0 if

(

x

,

y

)

∈ ∂

+

,

 >

0 if

(

x

,

y

)

∈ 

0

\

+

,

(21)

Fig. 4. Thelevelsetfunction(coloured)iszeroattheblackline,∂+.Theglacierboundaryconsistsofthepartsofthisblacklineandthemountain (dashedline)whichliesinbetweentheirintersection.

and

|∇|

=

1 at

∂

₊.Thedomain



₊ isadomainwhichislargerthantheglacierdomain



,andwhoseupperboundary coincides with



s,seethe blacklineinFig. 4.As theglacierdeforms, thelevel setfunction isconvectedby thevelocity field accordingto thefree surface equation (4). Toaccount foraccumulation/ablationtermas,we usethat as

:= −

a⊥s

∇

andsimplyaddanextracomponentuacc

= [

0

,

as⊥

]

T tothevelocitysothat(4) canbewrittenas

∂



∂

t

+ (

u

+

uacc

)

· ∇ =

0 in



0

.

(22)

Anewglacierboundary

∂

isthendefinedinthefollowingway.Theintersectionisfoundbetweenthelinewhere



=

0 (the black lineinFig. 4) andthe underlyingmountain(the dashed line inFig. 4). Anew



s isthen definedasthe part of theline



=

0 that liesin betweenthe intersectionpoints. In the sameway a new



b isdefined asthe partof the mountainboundarythatliesinbetweentheintersectionpoints.Inthiswaythemargins(thefrontandback)oftheglacier isfreetomove.The reasonweuseagreatdomain



₊ insteadofsimplyhavingalevelsetdefinedbyonlytheextension of



sispurelytechnical,i.e.itisbecausethecodeiswritteninsuchawaythattheline



=

0 isassumedtobe aclosed polygon. Inpractice the mountainboundary wouldbe given by databut inthis study we find a curve representingthe mountain, weextrapolate



b sothatit isprolongedtothesides. Ifthebaseoftheicewouldalsobe moving,e.g. dueto subglacialchannels,itispossibletorepresent



b withalevelsetfunction.Analternativeapproachistodefinethewhole glacier boundaryusing a levelset function, by replacing



₊ by



in(21), butwe found doing so smoothens the sharp edgesslightlyatthefrontandback,leadingtoartificialmassloss.

Note that asthe full Stokes equations are only solved in



h, we need toextend the velocity to



0. Thisis done by solvingtheLaplaceequationinthearea outsidetheglacierwiththecomputedvelocityasaDirichletconditionat

∂

and ahomogeneous Neumannconditionat

∂

0 (i.e.atthedashed lineinFig.2).Inordertosettheboundaryconditions,the CutFEMmethodisagainused,withanactivedomainconsistingofthecutelementstogetherwithelementscompletely out-sidetheglacier(thegreycolouredelementsinFig.2).Thecomputationalworkforsolvingthisproblemissmallcompared tosolvingthefullStokesequations.Thecompletevelocityfieldforallof



0isconstructedbycombiningthetwosolutions. Theextendedvelocityinnodeswithintheactiveglaciermesh



h

(

t

)

(i.e.inthemagentaandcyanareaofFig.2)aregiven bythefullStokessolutionandthesolutioninthenodesoutsidetheglacier(i.e.belongingonlytothegreycolouredareain Fig.2)aregivenbythesolutiontotheLaplaceequation.

In ordertoretain an accurate representationofthe geometry,equation (22) is solved ona meshtwice asfineasthe meshusedtosolve thefullStokesequations,employing linearinterpolation toobtainthevelocityuh/2 onthefinemesh. Todiscretise(22) weuseCrank-Nicolsonintimeandstabilisedlinearfiniteelementsinspace. Thediscreteproblemreads: Findthelevel-setfunction



h/2

∈

X0,h/2 suchthat



n_h+_/21

− 

n_h/2



t

, ξ

h/2

+

1 2

(

u

ˆ

n h/2

· ∇(

nh+/21

+ 

n h/2

), ξ

h/2

)

(23)

+

1 2s

(

n+1 h/2

+ 

n h/2

, ξ

h/2

)

=

0

∀ξ

h/2

∈

X0,h/2

,

whereu

ˆ

=

u

+

uacc andthestabilisation termsisintroducedtocontroltransportinstabilities

s

(

h

, ξ

h

)

= β





F∈Fh()

h2_n

|ˆ

uh

·

nF

|



nF

· ∇

h



,



nF

· ∇ξ

h





F

,

(24)

where

β

isauserdefinedparameter.Thisstabilisation termisthesameasin[22],exceptforthatweusethescalinghn insteadofh sincethisprovedmorestableinexperimentsonanisotropicmeshes.

Inordertoidentifythenew



s,i.e.thelinewherethelevelsetfunction



iszeroandthuschangingsign,itisimportant that

|∇|

remains closeto1.Therefore theadvectionof



is followedbya reinitialisation step.Availablereinitialisation

Table 1

Parametervaluesusedinthenumericalexperiments.

Description Notation Unit square Arolla glacier

Rate factor (MPa−3_year−1_{) A 2}1−1p ₁₀₀

Density (MPa/(m2_year2_{) ρ – 9}

.1380·10−19

Gravity (m/year2_{) g – [0}

,−9.7692·1015_]

Slip length  0 0,10300

Non-linear parameter p 4/3 4/3

Critical strain rate crit 10−7 ₁₀−7

Basal melt rate (m/year) a⊥_b 0 0, 0, 0.01

Surface accumulation (m/year) a⊥s 0 0, 0.5, 0

Tolerance, non-linear iter. tol 10−4 ₁₀−4

Mesh resolution [hx,hy] [2−i,2−i],i=2, . . . ,6 c 2−4[ √

8,√1

8],c=1,0.75

Inf-sup stability parameter α 0.1 10−6

Ghost penalty, velocity βu 0.1 1.0

Ghost penalty, pressure βp 0.1 0.01

Transport stabilisation par. β - 0.01

Tangential Nitsche penalty γt 0.05 0.005

Normal Nitsche penalty γn 0.05 0.005

techniquestendtosmoothverysharpedges,butthisiscircumventedbyintroducing



+.Toensurehighaccuracywekeep thefinemesh

T

0,h/2 andemploythelocalprojectionreinitialisationstrategy of[23] attheelements whicharecutby the boundary. The localprojection reinitialisation isbased onfinding an exact lineardistancefunction ineach element. This proceduregivesagloballydiscontinuousspacewhichisthenprojectedontoacontinuouspiecewiselinearfunctions.This methodworks alsoforunstructuredmeshes. Atelements that are notcut by theboundary thefastsweeping methodof [24] isused toextendthisdistancefunctioniteratively.The finalnewboundary is

∂

usedto define anewactive mesh

T

h

()

inthenexttimestep.

5. Numericalexperiments

We conduct fournumerical experiments,the purposeof whichis toshow that CutFEM can be usedto solve thefull Stokesequationsandthatitisausefulmethodforsimulationsoficedynamics.

Experiment 1a A unit square test inorder to verifythat thesolution converges atthe expected rateforthe full Stokes equationsandisrobustwithrespecttopenaltyparameters.

Experiment 2a AsteadysimulationoftheArollaglacier,withnoslipconditionsatthebase.

Experiment 2b AsteadysimulationoftheArollaglacierwithpartialslipconditionsattheicebed,representingasubglacial lakefoundundertheglacier.

Experiment 2c AtimedependentsimulationofArollaglacier,simulatingthefreesurfaceevolving betweenyear1930and 1932.

ThevaluesofallnumericalandphysicalparametersusedintheseexperimentsaregiveninTable1.

5.1. Experiment1- convergenceandrobustnessonaunitsquare

WefollowExample2insection4.8in[25].Thedomainisaunitsquare



:= (−

0

.

5

,

0

.

5

)

× (−

0

.

5

,

0

.

5

)

andtheanalytical solutionisgivenby ux

= (

x2

+

y2

)

a−1 2 _y

_,

₍₂₅₎ uy

= − (

x2

+

y2

)

a−1 2 _x

_,

₍₂₆₎ p

= (

x2

+

y2

)

b2_xy

−

_C

_.

₍₂₇₎

Therighthandsideofequation(1) isf

= ∇ ·

S

(

Du

)

+ ∇

p insteadof

ρ

g andDirichletconditionsuD

=

u

|

∂ areprescribed on the entire boundary

∂

, i.e.



=

0 in (7). To obtain a unique pressure we also enforce the condition



_pd



=

0, implying that C

=



_

(

x2

+

y2

)

b2_xyd

_

_{in (25). The ratefactor A issetso that} 1

2A

1−p

₌

_{1 inorder tobe consistent with} [25]. The constants a

,

b

∈ R

determine the regularity of the solution. The solution should be regular enough to satisfy

F(

Du

)

:= (



+ |

Du

|)

p−22 _Du

∈

_W1,2

_()

2×2 _andp

∈

_W1,p_{,whichcorrespondstoa}

_>

_{1 andb}

_>

−

₂

_/

_p

−

_{1 [16,17,25].Inorder} totestthemethodunderasdifficultconditionsaspossiblewechosea

=

1

.

01 andb

= −

2

/

p

−

0

.

99,asin[25].Giventhese regularityassumptions,optimalconvergenceestimates

Fig. 5. Results of Experiment 1 - the unit square test.

||

u

−

uh

||

1,p

≤

cuh

,

||

p

−

ph

||

p

≤

cph2/p



,

(28)

hasbeenestablishedinthecaseof

p

∈ (

1

,

2

)

andinf-supstableelements orfirstorderequalelements stabilisedwiththe localprojectionstabilisation,forsomeconstantscu andcp [16,17,25].InthisexperimentwetestnumericallyiftheCutFEM solutionappeartofulfilsimilarconvergenceestimatesandhowsensitiveitistothechoiceofthepenaltyparameters

γ

D,

β

u,and

β

p.

Toensurethatcellsarecut,thebackgroundmeshisrotatedwithrespecttothesquarewithanangleof

√

2

π

/

4 radians. For the convergence experiment we then run simulations on a series with meshes with a cellsize hx

=

hy

=

2−i

,

i

=

2

,

3

,

4

,

5

,

6.TheresultisshowninFig.5b.The velocityconvergesasexpectedandthepressureconvergesfasterthanthe expectedrateofh0.5_.

Forthesensitivitystudy wesethx

=

hy

=

2−3 andrunthreesetsofexperiments.Inthefirstwekeep

β

u and

β

p asin Table1andlet

γ

D takethevalues10−1

,

1

,

10

,

102

,

103,and104.Inthesecond wekeep

γ

D and

β

u asinTable1andlet

β

p take thevalues10−6

,

10−3

,

10−1

,

102

,

106,and 108 andinthethird wekeep

γ

D and

β

p asinTable1andlet

β

u take the values10−10

,

10−6

,

10−1

,

102

,

1010

,

1011

,

5

·

1010,and 1012.The resultsare showninFig.5c.The accuracyis notvery sensitivetochoiceofpenaltyparameters.

5.2. Experiment2a- theArollaglacierwithnoslipconditions

Inthisexperimentwecomputethevelocityandpressurealonga5 kmlongcentralflowlineoftheHautGlacierd’Arolla, asitwas measuredduring theLittle IceAgein1930(see Fig.1).The surfaceandbedrocktopography isavailableonthe ISMIP-HOM website[26] witharesolutionof100 m.Thephysicalparameters

ρ

,g and A aresetaccordingtothe ISMIP-HOMbenchmark,see Table1.Thebackgroundmesh consistsofrectangles ofwidthhx

=

2−4

√

8 andheight hy

=

2−4 1√₈, andwescaletheproblemsothatthelengthunitiskminsteadofm,inordertoensurethathx

<

1.

ThecomputedvelocitycomponentsandpressureareseeninFig.6.TheresultsareinagreementwithotherfullStokes modelsintheISMIP-HOMbenchmark,see[15].

Fig. 6. Results of Experiment 2a - Arolla with no slip conditions.

5.3. Experiment2b- theArollaglacierwithpartialslipconditions

Inthisexperiment,wechangethebasalboundaryconditionssothat





=

10300 _{if 2200 m}

_≤

_x

_≤

_{2500 m}



=

0 elsewhere (29)

TheseboundaryconditionsrepresentssubglacialwaterfoundundertheArollaglacierasdescribedintheISMIP-HOM bench-mark.ThelakeismarkedblueinFig.1.Thenumber10300_{ischosenhighenoughtoeffectivelymodel}

_

_{→ ∞}

_,sothatthere is free slipbetween the iceand the underlying water. The boundary conditions teststhe implementation of the Navier slipconditionsandarechallenging astheyare discontinuous.In ordertohavemorethanone meshnode residinginthe narrowregionoffreeslipwehererefinethemeshslightlyincomparisontoExperiment2a,sothathx

=

0

.

75

·

2−4

√

8 and hy

=

0

.

75

·

2−4 1√₈.

TheresultsshowninFig.7areinagreementwithotherfullStokesmodelsintheISMIP-HOMbenchmark,see[15] and [27], giventherelatively coarse resolutionhere. The coarseresolution usedinthispaperwas chosen to demonstratethe abilitytorepresentgeometrywithoutafinemesh.Inpracticeonecanrefinethebackgroundmeshinselectedareas.

5.4. Experiment2c- timeevolutionofthe Arollaglacier

Inthe finalexperimentwe let theicesurface evolvefromyear1930toyear1932. Thisexperimentis notpartofthe ISMIP-HOM benchmark,andwe insteadcomparetotheresultsin[28]. Weapply noslipconditionsatthebaseasthisis thecasein[28].The timestepis0.2years.Thesurface aftertwoyears(in1932) ispresentedinFig.8together withthe initialgeometryat1930.Theicethicknessdecreasesintheupperpartoftheglacierandtheiceaccumulatesinthelower part,i.e.theglacierismovingdownthemountain.Theresultsareinagreementwiththeresultsof[28].

In climate relatedapplications, it is important formodels to preserve mass. We therefore measure the area changes throughout the simulation, see Fig. 8b. Some area is artificially lost (22

.

3 m2_{), corresponding to about 2 mm artificial} surfacemelt peryear,oranrelativeareachange

|

A

(

t

=

2

)

−

A

(

t

=

0

)

|/

A

(

t

=

0

)

ontheorderof10−5_{.Giventhefairlyhigh}

Fig. 7. Results of Experiment 2b - the Arolla glacier with partial slip conditions.

non-lineartolerance



tol,coarsemeshandthestabilised equalorderlinearelements,thisiswithin theexpectedrange.It canbecomparedtoarelativevolumechangeof10−6_{fortheicesheetmodelElmer/Ice,foranexperimentwithanicewith} constantviscosityisflowingoveranartificialsinusoidal bedrockuntilsteadystateisreached[29].

Finally,todemonstratethatis possibletomodelaccumulationas attheicesurface,aswellasbasalmelta_b⊥atthebase, weruntwomoresimulationsandmeasureareachange.Inthefirstsimulation,wemodelhalfameterofaccumulationper yearovertheentiresurfacebysettinga⊥_s

=

0

.

5.Sincethelengthoftheglacierchangesthroughoutthesimulationsslightly (asthefrontandbackmoves)itisdifficulttocalculateexactlyhowlargetheareachangeshouldbe.Onaveragetheglacier lengthishoweverabout4890

.

1 m long,meaningthattheglaciershouldintwoyearsgain4890

.

1 m

·

0

.

5 m/year

·

2 years

=

4890

.

1 m2.Thisagreeswiththemeasuredarea gainof4877

.

6 m2,consideringthatthenumericalerrorisabout20 m2 (see Fig.8b).It shouldbe notedthatin realitythe accumulationontheArollaglacieris about

−

2 m[30]. Inthesecond simulation,wesetabasalmelta⊥_b

=

0

.

01 m/year,i.e.wesetgd

= (

0

,

−

0

.

01

)

.Thisshouldresultinaarea loss of4890

.

1 m

·

0

.

01 m/year

·

2 years

=

97.8m2,which alsoagreeswell withthemeasured arealossof120

.

2 m2,againconsidering the numericalerror.

6. Summaryandconclusion

The unfitted sharp interface finiteelement method CutFEMhas been appliedto simulate glaciersin two dimensions. A strategy forhandling non-linear viscosities, thindomains, variableslip conditions andan evolving domain withsharp edges has beendeveloped.The order ofconvergenceis atleastashigh asexpectedfromfinite elementtheory for non-Newtonian fluids.The methodproduces accuratevelocity andpressureprofiles fortheSwissglacierHautGlacierd’Arolla, forfullyfrozenbasalconditionsaswellasforanunderlyingsubglaciallake.Thedeformationoftheglacierfromyear1930 to 1932has been simulatedwithand withoutbasal melt and surface accumulation,showing that the movement ofthe glaciersurfaceiswellmodelleddespitenotbeingrepresentedbymeshnodes.

The CutFEM method appears to be a viable option for ice dynamics simulations which avoids issues with meshing andremeshing. The CutFEMmethodhas beenimplemented inthreedimensions inmaturesoftwares such asdeal.II and FEniCS[31,32] sothatenextensiontolargescalethreedimensionalicesheetsimulationsispossiblebyfollowingthesame

Fig. 8. Results of Experiment 2c - transient simulations of the Arolla glacier.

procedure asinthispaper.Other interesting extensionsofthisstudy aremarine terminating glaciers,configurationswith topologicalchanges,aswellasperformancecomparisonstoclassicalfittedFEMonlargescaleproblems.

CRediTauthorshipcontributionstatement

JosefinAhlkrona: Conceptualization,Formalanalysis,Investigation,Methodology,Software,Visualization,Writing– orig-inaldraft. DanielElfverson: Conceptualization,Formal analysis, Investigation, Methodology, Software,Writing –review & editing.

Declarationofcompetinginterest

The authors declare that they haveno known competingfinancial interests or personal relationships that could have appearedtoinfluencetheworkreportedinthispaper.

Acknowledgements

TheauthorsacknowledgethefundingprovidedbytheSwedishe-ScienceResearchCentre(SeRC).Partoftheresearchwas conductedduringaresearchvisitatKielUniversity,fundedbytheClusterofExcellence80“TheFutureOcean”.The“Future Ocean” is funded within the framework of the Excellence Initiative by the Deutsche Forschungsgemeinschaft (DFG) on behalfoftheGermanfederalandstategovernments.WearegratefulforfruitfuldiscussionswithChristianHelanowandto theanonymousreviewerwhothoroughlyreadthemanuscriptandwhosesuggestionsofadditionalnumericalexperiments improvedthestudysignificantly.

References

[1]J.Church,P.Clark,A.Cazenave,J.Gregory,S.Jevrejeva,A.Levermann,M.Merrifield,G.Milne,R.Nerem,P.Nunn,A.Payne,W.Pfeffer,D.Stammer,A. Unnikrishnan,ClimateChange2013:ThePhysicalScienceBasis.ContributionofWorkingGroupItotheFifthAssessmentReportofthe Intergovern-mentalPanelonClimateChange,CambridgeUniversityPress,Cambridge,UnitedKingdomandNewYork,NY,USA,2013,pp. 1137–1216.

[2] H.-O.Pörtner,D.Roberts,V.Masson-Delmotte,P.Zhai,M.Tignor,E.Poloczanska,K.Mintenbeck,A.Alegría,M.Nicolai,A.Okem,J.Petzold,B.Rama, N.e.Weyer,IPCCSpecialReportontheOceanandCryosphereinaChangingClimate,2019,inpress,https://www.ipcc.ch/srocc/.

[3] E.Larour,H.Seroussi,M.Morlighem,E.Rignot,Continentalscale,highorder,highspatialresolution,icesheetmodelingusingtheIceSheetSystem Model(ISSM),J.Geophys.Res.,EarthSurf.117(2012),https://doi.org/10.1029/2011JF002140.

[4] N.Petra,H.Zhu,G.Stadler,T.J.R.Hughes,O.Ghattas,AninexactGauss-Newtonmethodforinversionofbasalslidingandrheologyparametersina nonlinearStokesicesheetmodel,J.Glaciol.58(2012)889–903,https://doi.org/10.3189/2012JoG11J182.

[5] D.J.Brinkerhoff,J.V.Johnson,Dataassimilationandprognosticwholeicesheetmodellingwiththevariationallyderived,higherorder,opensource,and fullyparallelicesheetmodelVarGlaS,Cryosphere7(2013)1161–1184,https://doi.org/10.5194/tc-7-1161-2013.

[6] O.Gagliardini,T.Zwinger,etal.,CapabilitiesandperformanceofElmer/Ice,anewgenerationice-sheetmodel,Geosci.ModelDev.6(2013)1299–1318,

https://doi.org/10.5194/gmd-6-1299-2013.

[7] E.Burman,S.Claus,P.Hansbo,M.G.Larson,A.Massing,Cutfem:discretizinggeometryandpartialdifferentialequations,Int.J.Numer.MethodsEng. 104(2015)472–501,https://doi.org/10.1002/nme.4823.

[8] T.-P.Fries,T.Belytschko,Theextended/generalizedfiniteelementmethod:anoverviewofthemethodanditsapplications,Int.J.Numer.MethodsEng. 84(2010)253–304,https://doi.org/10.1002/nme.2914.

[9] J.Nitsche,ÜbereinVariationsprinzipzurLösungvonDirichlet-ProblemenbeiVerwendungvonTeilräumen,diekeinenRandbedingungenunterworfen sind,Abh.Math.Semin.Univ.Hamb.36(1971)9–15,https://doi.org/10.1007/bf02995904.

[10] A.Hansbo,P.Hansbo,Anunfittedfiniteelementmethod,basedonNitsche’smethod,forellipticinterfaceproblems,Comput.MethodsAppl.Mech. Eng.191(2002)5537–5552,https://doi.org/10.1016/S0045-7825(02)00524-8.

[11] A.Massing,M.G.Larson,A.Logg,M.E.Rognes,AstabilizedNitschefictitiousdomainmethodfortheStokesproblem,J.Sci.Comput.61(2014)604–628,

https://doi.org/10.1007/s10915-014-9838-9.

[12] T.Frachon,S.Zahedi,Acutfiniteelementmethodforincompressibletwo-phaseNavier–Stokesflows,J.Comput.Phys.384(2019)77–98,https:// doi.org/10.1016/j.jcp.2019.01.028.

[13] G.Cheng,P.Lötstedt,L.vonSydow,AfullStokessubgridschemeintwodimensionsforsimulationofgroundinglinemigrationinicesheetsusing Elmer/ICE(v8.3),Geosci.ModelDev.13(2020)2245–2258,https://doi.org/10.5194/gmd-13-2245-2020.

[14] J.H.Bondzio,H.Seroussi,M.Morlighem,T.Kleiner,M.Rückamp,A.Humbert,E.Y.Larour,Modellingcalvingfrontdynamicsusingalevel-setmethod: applicationtoJakobshavnIsbræ,WestGreenland,Cryosphere10(2016)497–510,https://doi.org/10.5194/tc-10-497-2016.

[15] F.Pattyn,L.Perichon,A.Aschwanden,B.Breuer,B.deSmedt,O.Gagliardini,G.H.Gudmundsson,R.Hindmarsh,A.Hubbard,J.V.Johnson,T.Kleiner,Y. Konovalov,C.Martin,A.J.Payne,D.Pollard,S.Price,M.Rückamp,F.Saito,O.Sou˘cek,S.Sugiyama,T.Zwinger,Benchmarkexperimentsforhigher-order andfull-Stokesicesheetmodels(ISMIP-HOM),Cryosphere2(2008)95–108,https://doi.org/10.5194/tc-2-95-2008.

[16] L.Belenki,L.C.Berselli,L.Diening,M.R ˚užiˇcka,Onthefiniteelementapproximationofp-Stokessystems,SIAMJ.Numer.Anal.50(2012)373–397,

https://doi.org/10.2307/41582741.

[17] A.Hirn,Approximationofthep-Stokesequationswithequal-orderfiniteelements,J.Math.FluidMech.15(2012)65–88,https://doi.org/10.1007/ s00021-012-0095-0.

[18] M.Winter,B.Schott,A.Massing,W.Wall,ANitschecutfiniteelementmethodfortheOseenproblemwithgeneralNavierboundaryconditions, Comput.MethodsAppl.Mech.Eng.330(2018)220–252,https://doi.org/10.1016/j.cma.2017.10.023.

[19] E.Burman,A.Ern,Continuousinteriorpenaltyhp-finiteelementmethodsforadvectionandadvection-diffusionequations,Math.Comput.76(2007) 1119–1140,https://doi.org/10.1090/S0025-5718-07-01951-5.

[20]D.Boffi,M.Fortin,F.Brezzi,MixedFiniteElementMethodsandApplications,SpringerSeriesinComputationalMathematics,Springer,Berlin, Heidel-berg,2013.

[21] S.Frei,Anedge-basedpressurestabilizationtechniqueforfiniteelementsonarbitrarilyanisotropicmeshes,Int.J.Numer.MethodsFluids89(2019) 407–429,https://doi.org/10.1002/fld.4701.

[22]E.Burman,M.A.Fernández,Continuousinteriorpenaltyfiniteelementmethodforthetransientconvection-diffusion-reactionequation,ResearchReport RR-6543,Inria,FrenchInstituteforResearchinComputerScienceandAutomation,2008.

[23]N.Parolini,E.Burman,Alocalprojectionreinitializationprocedureforthelevelsetequationonunstructuredgrids,TechnicalReport CMCS-REPORT-2007-004,ÉcolePolytechniqueFédéraledeLausanne,2007.

[24] Y.-H.R.Tsai,L.-T.Cheng,S.Osher,H.-K.Zhao,FastsweepingalgorithmsforaclassofHamilton—Jacobiequations,SIAMJ.Numer.Anal.42(2003) 673–694,https://doi.org/10.1137/S0036142901396533.

[25]A.Hirn,Finiteelementapproximationofproblemsinnon-Newtonianfluidmechanics,Ph.D.thesis,Naturwissenschaftlich-Mathematischen Gesamt-fakultätderRuprecht-Karls-UniversitätHeidelbergundderMathematisch-PhysikalischenFakultätderKarls-UniversitätPrag,2011.

[26] F.Pattyn,ISMIP-HOMbenchmarkwebsite,http://homepages.ulb.ac.be/~fpattyn/ismip/.(Accessed 30June2019).

[27]C.Helanow,Effectsofnumericalimplementationsoftheimpenetrabilityconditiononnon-linearStokesflow:applicationstoicedynamics,e-prints, arXiv:1907.12315,2019.

[28] J.Ahlkrona,V.Shcherbakov,Ameshfreeapproachtonon-Newtonianfreesurfaceiceflow:applicationtotheHautGlacierd’Arolla,J.Comput.Phys. 330(2017)633–649,https://doi.org/10.1016/j.jcp.2016.10.045.

[29] O.Gagliardini,T.Zwinger,TheISMIP-HOMbenchmarkexperimentsperformedusingthefinite-elementcodeElmer,Cryosphere2(2008)67–76,https:// doi.org/10.5194/tc-2-67-2008.

[30] R.Dadic,J.G.Corripio,P.Burlando,Mass-balanceestimatesforHautGlacierd’Arolla,Switzerland,from2000to2006usingdemsanddistributed mass-balancemodeling,Ann.Glaciol.49(2008)22–26,https://doi.org/10.3189/172756408787814816.

[31] C. Gürkan,S.Sticko, A.Massing,StabilizedcutdiscontinuousGalerkin methodsfor advection-reactionproblems,SIAMJ.Sci.Comput. 42(2020) A2620–A2654,https://doi.org/10.1137/18M1206461.

[32] A.Massing,B.Schott,W.Wall,AstabilizedNitschecutfiniteelementmethodfortheOseenproblem,Comput.MethodsAppl.Mech.Eng.328(2018) 262–300,https://doi.org/10.1016/j.cma.2017.09.003.