Regions of high velocities after inversion

analuzardi

Hi,

I'm having problems with my inversion for basal friction coefficient for a Greenland-wide simulation. My initial guess for basal friction coefficient seems good:

My cost function J:

Here is md.results.StressbalanceSolution.Vel - md.inversion.vel_obs:

The calculated velocity seems good at most parts. But at some very specific locations (near ice margin, but in regions where there is ice), there are some values that are too high.

Difference between calculated friction coefficient and my initial guess (it didn't change much):

I already tried lots of things to solve this, like: increasing resolution, smoothing surface elevation and thickness, playing with values for penalizing for velocity difference, log velocity difference and gradients. Unfortunately, nothing seems to work.

I also plotted several variables to see if I could find something crazy at those regions with high velocity (e.g., thickness, geometry, bedrock, ice surface slope, pressure, etc.), but didn't find anything that seem to be wrong. I'm stuck 🙁

I would appreciate any tips of what could be going wrong or things to check.

Also, for a good inversion, how much should I get for RMSE between md.inversion.vel_obs and md.results.StressbalanceSolution.Vel?

Thanks,
-Ana

mathieumorlighem

Hi Ana

Good to see you on the forum!
Inversions are hard, and there is no recipe. Here is what I would recommend to do:

run a stress balance without inversion and plot the speed. Does it look good? or do you have some points with a very high speed? If you do, you will need to take care of that first (by changing the outline, mask, etc).

if you have a good initial speed, you should try to run the inversion for just ₁₀ steps and look at the gradient, which tells you where the inversion is trying to change the friction:

plotmodel(md,'data', md.results.StressbalanceSolution.Gradient1);

Hopefully these two steps will help you get started
Cheers
Mathieu

analuzardi

Hi Mathieu,

Thanks for the prompt response!

I did what you suggested, I ran the stress balance with md.inversion.iscontrol=0. The velocity seems ok at most parts:

But I have crazy values at some regions:

I'll try changing the outline and mask tomorrow and see how it goes. I'll keep you posted.

Thanks,
-Ana

mathieumorlighem

yeah, these will definitely be problematic... good luck!

analuzardi

I have an update! After spending some time on it, I realized my ice mask is messed up. After interpolating the ice mask from grid to mesh, I get something like this:

I was careful to prepare my ice mask (that comes from BedMachine) before giving to ISSM. I deleted surrounding islands and got rid of the holes in the ice sheet interior. However, the above messy ice mask is what I get after this line of code:

mask_mesh = InterpFromGridToMesh(x_mask.astype(float), np.flipud(y_mask.astype(float)), np.flipud(mask.astype(float)), md.mesh.x, md.mesh.y, 0)

and I'm using the output of this to feed md.mask.ice_levelset. So I think I need to fix this first. That's probably what's causing issues with my inversion. I tried setting the mask with "setmask" function with an exp file so that I wouldn't have to do this interpolation, but that function didn't work for me. Would you have any tips on how to fix the mask after interpolating? Or if it's best to define the mask from an exp file instead?

Thanks,
-Ana

mathieumorlighem

Hi Ana

have you tried something like that (to translate into python):

mask = interpBedmachineGreenland(md.mesh.x,md.mesh.y,'mask');
md.mask.ice_levelset= -1*ones(md.mesh.numberofvertices,1);
pos = find(mask<1);
md.mask.ice_levelset(pos)=1;

otherwise... I still think you would save yourself months (if not years!) of work by starting from Josh's model.
My 2 cents
Mathieu

analuzardi

Hi Mathieu,

The function interpBedmachineGreenland is not installed in the EIS SMCE pCluster 🙁 I used to be able to find some ISSM Python functions in this link:

http://icewiki.umt.edu/index.php/ISSM_Python_assimilation

and then upload it to EIS SMCE pCluster. But the link is not working anymore. I also tried finding this function on my local ISSM installations but couldn't find it. Would you know how I can get this function (hopefully in python)?

PS: I really want to try to have my own set up since I would then understand every line of it and that would be useful for playing with initialization, which is the point of my research. Hopefully with a lot of effort and help I can make it work, although I know it takes more time. But Iet's se.. I can still change my mind if I get really really stuck. I have hope!

Thanks,
-Ana

mathieumorlighem

Hi Ana

I worked with NSIDC a few years ago and there is a python equivalent here: https://github.com/nsidc/nsidc0756-scripts/blob/main/interp_bedmachine_antarctica.py (you just need to change Antarctica to Greenland, but the rest of the function should work out of the box)
I hope this helps!
Mathieu

analuzardi

Awesome thanks!! I'll try using this function and then keep you posted here whether it works.

-Ana