Hi folks:

Is it possible to transiently change the upperwater_melt_rate and deepwater_melt_rate parameters when using the linearbasalforcings() function? For example, if I wanted to have a deepwater_melt_rate of 20 m/yr for the first 10 years, and then change to 30 m/yr for the next 10 years, is this possible?

Does this do what I hope?
md.timestepping.start_time = 0;
md.timestepping.final_time = 30;
md.timestepping.time_step = 1/12;
md.basalforcings.deepwater_melt_rate = [20, 30; 0, 10+(1/12)]

I think this approach would linearly interpolate between 20 m/yr at time 0 and 30 m/yr at time 10+(1/12)? If I want a step-change to occur, do I just need to explicitly define a value for each time step (at leat between time 0 and 10+(1/12)?

Cheers,

Lawrence

Hi @lawrence-bird

Good news! It is already implemented and should work just as you described. If you want a step change you have two options. Either you add one more "time slice" to your time series:

md.basalforcings.deepwater_melt_rate = [20, 20, 30; 0, 10, 10+(1/12)];

so that the transition happens quickly between 10 and 10+1/12, or you can set md.timestepping.interp_forcing = 0 and none of the time varying inputs or parameters in the model would be interpolated in time.
Best
Mathieu

    a month later

    Hi @mathieumorlighem

    I'm hoping to clear something up RE: the way that ISSM applies forcing during transient runs.

    In this example, are the basal melt rates applied to the current timestep, or the next timestep? For example, are the effects of the melt rate of 30 m/yr seen at timestep 10+1/12, or 10+2/12? I am looking to extract the results at a given timestep when my basal melt rate changes and want to make sure that I'm getting the correct timestep.

    Similarly, does ISSM apply the md.levelset.spclevelset in the same way? I'm applying a change to the md.levelset.spclevelset at time 10, but I don't see any changes to the velocity field until 10+(1/12).

    I interpret this that ISSM uses the md.levelset.spclevelset at time 10 to compute a forward step and the results of this change are returned at 10+(1/12). Is this true? Is basal melt applied in the same way?

    Thanks!

    Lawrence

    Hi @lawrence-bird

    in general, all transient fields are "linearly" interpolated in time. For example, if you specify the melt to be 10 m/yr at t=0 yr, and 30 m/yr at t= 1 yr, the model will use a melt of 15 m/yr for time t = 0.5 yr. If you want to visualize it, you should add BasalforcingsFloatingiceMeltingRate to your list of requested output.
    If you do not want a linear interpolation but a step change instead, you need to change md.timestepping.interp_forcing to 0. Same thing for spclevelset (which should be defined as a distance field, not ±1!).
    I hope this helps!
    Mathieu

      Hi mathieumorlighem

      Yep, this makes sense - thanks!

      But, to confirm, are the effects of forcings used at Time = t seen at Time = t+1? Does ISSM use the forcings at Time = t to fun one step forward, such that the results of these forcings are stored in t+1?

      As for the spclevelset, what do you mean when you say "should be defined as a distance field"? I ended up explicitly defining the levelset at each timestep because I couldn't get the interpolation to work as I wanted. But that's probably because I was using +/- 1 masks. How would you recommend going about this?

      Cheers,

      Lawrence

        The effects prescribed at time t are applied at time t.

        For the signed distance you can use reinitializelevelset. It will preserve the 0 contour but will make your levelset a distance field. You can convert any level set like this:

        for i=1:size(md.levelset.spclevelset,2)
   %Extract levelset for time i
   levelset = md.levelset.spclevelset(1:end-1, i);
   %Make signed distance
   levelset=reinitializelevelset(md, levelset);
   %Insert back into series
   md.levelset.spclevelset(1:end-1, i) = levelset;
end

        Mathieu

          a month later

          Hi mathieumorlighem:

          Thanks for the tip on the reinitializelevelset() option -- I'll have a look at this.

          I think I did a terrible job explaining my query about the timestepping! What I wanted to confirm is that the results of effects prescribed at Time t exist in time t+1 in the results? For example, if I have a timestep of 1-year (for simplicity), and I prescribe md.levelset.spclevet at time = 1, this levelset shows up in md.results.TransientSolution.MaskIceLevelSet at time =2. Similarly, the a basal melt rate prescribed using linearbasalforcings() at time = 4 will show up in md.results.TransientSolution.BasalforcingsFloatingiceMeltingRate at time = 5. Because one forward timestep is completed with the forcing at time t to generate the results at time t+1.

          Cheers,

          Lawrence

            that should not be the case... here is an example:

            md.smb.mass_balance = repmat(  [1. 2. 3. 4.],[md.mesh.numberofvertices, 1]);
md.smb.mass_balance(end+1,:) = [1. 2. 3. 4.];
md.timestepping.start_time = 0;
md.timestepping.final_time = 5;
md.timestepping.time_step  = 0.5;

md.timestepping.interp_forcing = 1;
md=solve(md,'Transient');

smb = cell2mat({md.results.TransientSolution(:).SmbMassBalance});
times = [md.results.TransientSolution(:).time];

[smb(1,:); times]

            will show this:

            SMB:    1.0000    1.0000    1.5000    2.0000    2.5000    3.0000    3.5000    4.0000    4.0000    4.0000
TIME:   0.5000    1.0000    1.5000    2.0000    2.5000    3.0000    3.5000    4.0000    4.5000    5.0000

            and if you use md.timestepping.interp_forcing = 0

             SMB:   1.0000    1.0000    1.0000    2.0000    2.0000    3.0000    3.0000    4.0000    4.0000    4.0000
 TIME:  0.5000    1.0000    1.5000    2.0000    2.5000    3.0000    3.5000    4.0000    4.5000    5.0000

            as expected
            Mathieu

            15 days later

            Hi @mathieumorlighem
            Thanks for this example. I played a bit with one of the nightly tests and adding in basal forcing and adjusting the spclevelset.

            I see that the md.results.TransientSolution.BasalforcingsFloatingiceMeltingRate and md.results.TransientSolution.MaskIceLevelset do appear to show up in the correct time slots (exactly as I provide them). Is it correct to assume that the effects of these changes on other fields (e.g., Velocity and Thickness) should also appear in the same time slot, or would there be a 1-timestep offset in these fields?

            Essentially, I'm trying to isolate the instantaneous effects of a calving event (and changes in basal melt) and I am trying to confirm whether I should be extracting the same timestep as the forcing field, or the next timestep from things like md.results.TransientSolution.Vel and md.results.TransientSolution.Thickness.

            I convinced myself that there was a 1-timestep offset here between the forcing fields and the resulting changes in velocity/thickness, but now I'm not so sure!

            Cheers,

            Lawrence

            Hi Lawrence

            this is a good point. The easiest is that you take a look at src/c/cores/transient_core.cpp and the function transient_step:

            at each time step i (time t):

            • compute new temperature
            • smb
            • hydrology
            • stress balance (speed)
            • move the ice front (according to the level-set)
            • new thickness based on new ice front and new speed
            • update grounding line
            • etc
            • save results for this time step

            what you can see is that the speed will use the PREVIOUS level-set (i-1) but the thickness update happens AFTER the level-set update (i). I hope this answers your question?
            Cheers
            Mathieu

              Thanks mathieumorlighem -- this helps!

              It looks like the masstransport (specifically, bmb_core -- I presume this is where basal melting is calculated?) occurs after the stress balance for a given timestep.

              So, similarly to the levelset changes, does this mean that the velocities calculated at (i) are based on the basal melt forcings specified at the previous step (i-1), but the thickness at (i) is based on the basal forcing at (i)?

              Cheers,

              Lawrence

                The ice velocity should not depend on the melt rate 🙂 the melt only affects masstransport_core.
                But, yes, the ice velocity at i is calculated as a function of the ice thickness and levelset from i-1
                Mathieu

                Of course. Okay, I think I'm putting this to bed! Knowing that there is an offset in the transient levelset (and not in the basal melt/mass transport) is helpful!

                Thanks @mathieumorlighem !

                Sorry @mathieumorlighem -- one last question on this!

                Given the timestep offsets, is there any way of specifying exact timesteps for the output frequency? I recognise that setting md.settings.output_frequency = 1 would give me all the timesteps I could ever want. But I'm wondering if I can define a series of output timesteps, something like (for example) [0.083, 0.5, 1, 1.083, 1.5, 2.0, 2.083] to reduce the file size but ensure I'm capturing all the timesteps I am interested in?

                Cheers.

                no unfortunately you can't (or you would need to hardcode it)... but you can change the number of requested_outputs to only save what you need to save