Greenland with SIA

andreasplach

Hello everyone,

I am trying to do a transient run for Greenland with SIA (for a start with SMB=0 everywhere) just to see how the various approximations behave. I am basically using the ISSM Greenland tutorial. I did transient simulations for the same domain with SSA and HO and they work fine, but when I try to do the same with SIA I get enormous peaks (>1e8m/yr) in the surface velocity after a couple of decades and the simulation stops after a while with:


 ??? Error using ==> ./solutionsequences/solutionsequence_linear.cpp
 solutionsequence_linear error message:    solver residue too high!: norm(KU-F)/norm(F)=0.000165554

Since the inversion seems not to be supported for SIA, I need to provide proper parameters for the friction law I use. I tried runs with the default friction law (md.friction=friction()) and the Weertman friction law (md.friction=frictionweertman()). All combinations I tried, gave similar results concerning the velocity peaks. Here are some of the friction combinations I tried:

1) Just the parameters from the Greenland tutorial:


md.friction=friction()
md.friction.coefficient=30*npy.ones((md.mesh.numberofvertices,1))
pos=npy.where(md.mask.groundedice_levelset<0)
md.friction.coefficient=0  # no friction applied on floating ice
md.friction.p=npy.ones((md.mesh.numberofelements,1))
md.friction.q=npy.ones((md.mesh.numberofelements,1))


md.friction=frictionweertman()
md.friction.C=1e-12*npy.ones((md.mesh.numberofvertices,1))
md.friction.m=1*npy.ones((md.mesh.numberofelements,1))

3) Parameter values I found in the Elmer/Ice online documentation.


md.friction=frictionweertman()
md.friction.C=2.412579e-2*npy.ones((md.mesh.numberofvertices,1))
md.friction.m=1.0/3.0*npy.ones((md.mesh.numberofelements,1))

For the ice rheology I just use the the values from the Greenland tutorial:


md.materials.rheology_n=3*npy.ones((md.mesh.numberofelements,1))
md.materials.rheology_B=paterson(md.initialization.temperature)
md.friction.q=npy.ones((md.mesh.numberofelements,1))
md.friction.p=npy.ones((md.mesh.numberofelements,1))

The whole domain is set to being grounded and ice is existing everywhere (minimum of 1m ice), only the edges are set to be the ice front:


pos0=npy.where(md.geometry.thickness<=0)
md.geometry.thickness=1
md.geometry.surface=md.geometry.thickness+md.geometry.base
md.mask.ice_levelset=0

But I also did a run where I did not set the ice front => same result. The velocity peaks appear to be at the edges of the domain or at least close to the edges.

In the end I would like to use a combination of SIA and SSA for Greenland, e.g. use SSA where surface velocity > threshold, else SIA. I hope that should be possible in theory?

Any ideas or recommendations are highly appreciated.

I am using the ISSM developer version at revision 20920 with the python interface.

Cheers

Andreas

mathieumorlighem

Hi Andreas,

welcome to the world of SIA modeling! Everything is local so if your slopes are locally high, the velocity might be very large and that will break the transient run, which is why we stay away from it.

SIA only includes vertical shear, so boundary conditions like water pressure at the ice front are not taken into account. Same for the friction law, we cannot impose the stress but the basal velocity only, whether you use weertmenfriction or friction you will get the same results since none of these fields are used by SIA.

Since you have the development version, you can open src/c/analyses/StressbalanceSIAAnalysis.cpp and look at StressbalanceSIAAnalysis:. You will see that several friction "laws" (which basically specify what the basal velocities should be) are implemented by the one that is activated by default is frictionlaw = 1:

case 1:


   /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
   ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope;

vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope;

You can try to change that friction law to 2 or 3 (and use md.friction=friction() and provide md.friction.coefficient that is used by these two laws). You can also try and smooth the surface elevation first so that your initial driving stress is no too high.

I hope that helps

Best

Mathieu

andreasplach

Hei Mathieu,

thank you for your fast reply! After playing around a bit with the coefficients, my simulation works with law 3. I am planning to use SSA for the fast flowing parts of my ice sheet in the future.

Are you aware of anybody using ISSM in a combination of SIA and SSA? There is this on entry about SIA/SSA in the forum, but there is no reply, so I do not know how that worked out.

https://issm.ess.uci.edu/forum/viewtopic.php?f=13&t=91

Cheers

Andreas

mathieumorlighem

Hi Andreas,

it is fairly easy, you could do that (for example):

md=setflowequation(md,'SSA',md.inversion.vel_obs>100,'fill','SIA');

That uses SSA in the regions where vel>100 m/yr and SIA everywhere else. Try it!

Mathieu

PS: we rarely use model coupling actually, so nobody is an "expert" in the field, you can also contact Helene Seroussi if you have further question about coupling.