Transient mask updates

lawrence-bird

Hi folks:

I'm wondering what the best approach is to implement periodic calving events within ISSM.

If I have defined ice front positions for given dates/events is it possible to adjust the md.mask at given timesteps within a transient model run to convert floating ice to ocean cells? Or would I have to run a transient simulation to the timestep prior to the calving event, save the model, adjust the mask, then reinitialize the model with the previous timestep and begin a new transient run?

I understand the ice front position can be adjusted using the modelextract function (as per Pig Sensitivity example), but I am not sure how best to implement this within a transient simulation.

Thanks for any tips!

Lawrence

mathieumorlighem

Hi Lawrence

the easiest is to create a time series for md.levelset.spclevelset and turn on md.transient.ismovingfront = 1. You will need to constrain every node in md.levelset.spclevelset (i.e. do not leave any NaN) and ISSM will interpolate in between time steps (for example, if you set a node as ice for t=0 and ocean for t=1, ISSM will figure out what to do for t=0.5). Here is an example for one of the nightly run tests (test208.m):

%Construct levelset time series (here we collapse the entire ice shelf at time t=1 yr)
md.levelset.spclevelset = [md.mask.ice_levelset md.mask.ice_levelset; 0 1];
pos = find(md.mesh.y>500e3);
md.levelset.spclevelset(pos,2) = +1;

%Set transient parameters and time stepping
md.transient.ismovingfront = 1;
%md.levelset.kill_icebergs = 0; %not needed for this experiment
md.timestepping.start_time = 0;
md.timestepping.final_time = 2;
md.timestepping.time_step  = .3;

%Set other parameters needed by movingfront_core (but won't be used because it is completely constrained)
md.calving.calvingrate         = zeros(md.mesh.numberofvertices,1);
md.frontalforcings.meltingrate = zeros(md.mesh.numberofvertices,1);

Let me/us know if you need help!
Best
Mathieu

lawrence-bird

Thanks mathieumorlighem!

I'll have a play with this and let you know how I go.

Is it possible to adjust/turn off the interpolation between timesteps? For example, if I wanted one large calving event at a give time step vs. a linear retreat of the ice front between prescribed ice front positions. Presumably the best approach for this would be to provide a complete md.levelset.spclevelset that explicitly includes what I'm looking for at each timestep (rather than letting ISSM do any interpolation between timesteps)?

Cheers,

Lawrence

mathieumorlighem

Hi @lawrence-bird

yes you can play with md.timestepping.interp_forcing (default is 1 = interpolate in between given times, 0 = do not interpolate). Note that this will apply to all forcings, SMB, melt etc...
Cheers
Mathieu

lawrence-bird

Hi @mathieumorlighem:

I'm just getting back to thinking about this a little more. How does the md.levelset.spclevelset object compare with the md.mask.ice_levelset and md.mask.ocean_levelset objects? Does the md.levelset.spclevelset only define ice/ocean, and md.mask.ice_levelset is used to define which of those ice nodes are grounded/floating?

Cheers,

Lawrence

mathieumorlighem

Hi Laurence

md.levelset.spclevelset only constrains the ice levelset (not the ocean levelset, which is controlled by hydrostatic equilibrium). If you run a model with movingfront on, then the level set equation needs to be solved to reset the ice level set based on the calving rate and ice speed. The trick you are using here is to constrain it everywhere with md.levelset.spclevelset instead of using a calving law. In other words, md.levelset.spclevelset is the time series that you want to apply to md.mask.ice_levelset throughout your simulation. I hope this makes sense?
Mathieu

lawrence-bird

Hi mathieumorlighem !

I hope it's okay to resurrect this old thread rather than make a new one!

I have been playing with this, but I'm running into some issues. I have a transient model that works fine when md.transient.ismovingfront = 0. However, when I set md.transient.ismovingfront = 1 and specify md.levelset.spclevelset for each timestep, the model fails on the second iteration.

I get this error and I'm not sure what this might mean, or how best to go about rectifying it...

call computational core:
iteration 1/300  time [yr]: 0.08 (time step: 0.08)
   computing new velocity
   computing slope of levelset function...
   computing calving and undercutting
   computing new ice front position
   extrapolation of Vx:
   extrapolation of Vy:
   extrapolation of Thickness:
   call computational core:
   looking for icebergs to kill
   reinitializing level set after killing 3 icebergs
   updating boundary conditions...
   computing basal mass balance
   computing mass transport
   call computational core
   computing new grounding line position
   updating vertices positions
   computing transient requested outputs
   saving temporary results
iteration 2/300  time [yr]: 0.17 (time step: 0.08)
   computing new velocity
solver residue too high!: norm(KU-F)/norm(F)=0.00543691 > 1e-06 (md.settings.solver_residue_threshold)
WARNING: Solver failed, Trying Recovery Mode
solver residue too high!: norm(KU-F)/norm(F)=0.00543691 > 1e-06 (md.settings.solver_residue_threshold)

[5] ??? Error using ==> ./modules/Solverx/Solverx.cpp:32
[5] Solverx error message: Recovery solver failed...


[18] ??? Error using ==> ./modules/Solverx/Solverx.cpp:32
[18] Solverx error message: Recovery solver failed...

Themd.transient.levelset object is simply an array of ice masks (one for each timestep), with -1 where ice exists and -1 for ocean nodes. md.levelset looks like this (default values used):

Level-set parameters:
         stabilization          : 1               -- 0: no, 1: artificial_diffusivity, 2: streamline upwinding
         spclevelset            : (53149x300)     -- Levelset constraints (NaN means no constraint)
         reinit_frequency       : 5               -- Amount of time steps after which the levelset function in re-initialized
         kill_icebergs          : 1               -- remove floating icebergs to prevent rigid body motions (1: true, 0: false)
         migration_max          : 1000000000000   -- maximum allowed migration rate (m/a)
         fe                     : 'P1'            -- Finite Element type: 'P1' (default), or 'P2'

Any tips would be appreciated!

Cheers,

Lawrence

mathieumorlighem

Hi @lawrence-bird,

the stress balance did not converge. It is not the end of the world, but you should take a look at the results of step 2 and see what's going on (you probably have a couple of nodes going too fast). Can you do md.settings.solver_residue_threshold = 0.01 temporarily, re-run your transient model to step 2, and plot the velocity? That will give you an idea.
Cheers
Mathieu

lawrence-bird

Thanks mathieumorlighem -- it was indeed a few nodes going too fast. A result of a spurious mask update I believe!

Cheers,

Lawrence