Implementation about SHAKTI Model

wudipao88

I'm currently working on implementing the SHAKTI model, specifically focusing on three examples from your paper published in 2018 (10.5194/gmd-11-2955-2018). In the paper, the example 1 is available in ISSM userguide, but the source code for examples 2 and 3 is not provided. The results I obtained for examples 2 and 3 using my own code consistently differ somewhat from those in your paper. So, I was wondering if you could provide the source code for examples 2 and 3. I'd like to identify where my implementation diverges from yours.

Best regard

aleahsommers

Hi-- we're happy to help, and glad to hear that you're interested in using SHAKTI. I am looking through old scripts to see if I can find the setup for the simulations in the 2018 paper (sorry that they weren't all archived back then for easy access). I'll get back to you soon -- in the meantime, if you'd like to share any more details about your results, maybe we can figure out what's going on.

Thanks,
Aleah Sommers

aleahsommers

@wudipao88 -- Here is a script to set up the SHAKTI example with 10 moulins. I just had to change a few lines of syntax to make it compatible with current ISSM (some things have changed since 2018). Let me know if you have questions!

runme-shakti-10moulins.m

4kB

shakti-10moulins.txt

1kB

shakti-10moulins.txt

107B

wudipao88

aleahsommers Thank you very much for your quick reply and code. Referring to the parabola equation corresponding to the ice surface in your code, I have successfully reproduced the effect of Example 2 (Channelization with multiple moulins). At the same time, the effect of Example 3 (Seasonal variation and distributed meltwater input) I ran is still a little different from the effect in your paper. I guess this is still because the parabola of the ice surface I constructed is different from yours, so can you provide the surface geometric information about Example 3?

All the best!

aleahsommers

wudipao88 Great! I'm glad to hear the geometry sorted out the simulations. Try this for the seasonal distributed example in the 2018 SHAKTI paper:

md.geometry.base = zeros(md.mesh.numberofvertices,1); % Flat bed
md.geometry.surface= 6*(sqrt(md.mesh.x+5e3) - sqrt(5e3)) + 550;

Let me know if that helps!

wudipao88

aleahsommers Thank you for your geometrical information. However, the results I get based on this geometry are still different from yours (especially the hydraulic head and effective pressure), and I don't know where the problem lies. This is my script for Example 3 (Seasonal variation and distributed meltwater input). Can you check where the problem is? If you don't have time, can you provide me with your script? Thank you very much!

Best wishes!

example-3.m

8kB

example-3.txt

113B

aleahsommers

wudipao88 Looking back at my old scripts, I think the exaggerated seasonal cycle was imposed after a spin-up with steady low input. And I think this is what we did for the forcing -- using a small nonzero englacial input during the "winter" before and after the big meltwater cycle.

md.timestepping.time_step=3600/md.constants.yts;
md.timestepping.final_time=365/365;
time=0:md.timestepping.time_step:md.timestepping.final_time;
md.hydrology.englacial_input = zeros(md.mesh.numberofvertices+1,numel(time));
md.hydrology.englacial_input(end,🙂=time;

% Seasonal input -- long winter, but smooth transitions
md.hydrology.englacial_input((1:end-1),🙂=3.17e-8360024365;
for tt=1:length(time)
if time(tt)>=0.4 && time(tt)<0.7
md.hydrology.englacial_input((1:end-1),tt)=-492.75.cos(2pi/0.3.(time(tt)-0.4))+492.75+1; % m/yr
elseif time(tt)>=1.4 && time(tt)<1.7
md.hydrology.englacial_input((1:end-1),tt)=-492.75.cos(2pi/0.3.*(time(tt)-1.4))+492.75+1; % m/yr
end
end

Also, it looks like we used md.materials.rheology_B(🙂= (2.5e-25)^-1/3, which is different that what is in your script.

See if that helps!

wudipao88

Thank you very much. Based on the "md.materials.rheology_B" you provided, my simulation results have improved and are closer to your results. However, there are still some gaps between the two results, especially in the effective pressure and gap height. The following figure is the variation curve of the minimum, maximum and mean gap height in my simulation for one year. Compared with your results, the maximum gap height is smaller. So I still don't know where the problem lies. I hope you can give me further help. Thank you very much!

In addition, the current application area of the SHAKTI model is mainly Greenland. I wonder if SHAKTI can be applied to Antarctica (including outlet glaciers and inland grounded glaciers)?

aleahsommers

I'm trying to figure out exactly what I did back in 2018 for the paper, and the scripts are a bit messy (sorry that it was not cleanly archived back then). I think that we did a short spin-up with uniform input of 1 m/yr everywhere, before beginning the seasonal forcing cycle. This is also described in the 2018 GMD paper on page 2965: "The model is first run with steady distributed input of 1 m a−1 in a spin-up stage with a time step of 1 h (steady state achieved in 4 days)." Did you do that part?

One other thing that comes to mind is that the current version of SHAKTI in ISSM imposes a limit on gap height to not exceed 1m (both in the source code and in the current binary code, I believe). That could be leading to some differences, since it looks like we reported a maximum gap height exceeding 1m during the extreme seasonal forcing.

In answer to your query about Antarctica -- yes, it would be great to apply SHAKTI in Antarctica. So far, we have mostly been working in Greenland and starting to do smaller mountain glaciers as well. I would be very happy to help if you'd like to get started with SHAKTI in Antarctica, and I think it could be really useful.

wudipao88

Yes, except for the seasonal forcing cycle between 0.4 and 0.7 years, the englacial input is set to a uniform input of 1 m/y. Below is my script for parameter settings and seasonal water input.

%======================================= Geometrical information
md.geometry.base = zeros(md.mesh.numberofvertices, 1);
md.geometry.surface = 6 * (sqrt(md.mesh.x + 5e3) - sqrt(5e3)) + 550;
md.geometry.thickness = md.geometry.surface - md.geometry.base;
%=======================================

%================== Material Information
md.materials.beta = 7.5e-8;
md.materials.heatcapacity = 4220;
md.materials.rho_ice = 910;
%==================

%====================================================== Basal sliding speed
md.initialization.vx = 10^-6 * md.constants.yts * ones(md.mesh.numberofvertices, 1);
md.initialization.vy = zeros(md.mesh.numberofvertices, 1);
%======================================================

%============================= Flow law parameters
md.materials.rheology_B(🙂 = (2.5e-25) ^ (-1/3);
md.materials.rheology_n = 3 .* ones(md.mesh.numberofelements, 1);
%=============================

%================================================== Temperature
md.initialization.temperature = (273 - 20) * ones(md.mesh.numberofvertices, 1);
%==================================================

%========================= Friction Coefficient
md.friction = frictionshakti(md.friction);
md.friction.coefficient = 20 .* ones(md.mesh.numberofvertices, 1);
%=========================

%================================================= Hydrological parameters
md.hydrology = hydrologyshakti();
md.hydrology.head = 0.5 * md.materials.rho_ice / md.materials.rho_freshwater * md.geometry.thickness + md.geometry.base;

random_vector = randn(md.mesh.numberofelements, 1);
mean_value = 0.01;
std_dev = 0.01 * 0.01;
perturbed_vector = mean_value + std_dev * random_vector;
md.hydrology.gap_height = perturbed_vector;

md.hydrology.bump_spacing = 2 * ones(md.mesh.numberofelements, 1);
md.hydrology.bump_height = 0.1 * ones(md.mesh.numberofelements, 1);
md.hydrology.englacial_input = 0 * ones(md.mesh.numberofvertices, 1);
md.hydrology.reynolds = 1000 * ones(md.mesh.numberofelements, 1);
%=================================================

%======================================== Water input
md.timestepping.time_step = 3600 / md.constants.yts;
md.timestepping.final_time = 365 / 365;
time = 0 : md.timestepping.time_step : md.timestepping.final_time;
time = time';

md.hydrology.englacial_input = zeros(md.mesh.numberofvertices + 1, numel(time));
md.hydrology.englacial_input(end, 🙂 = time;

% Seasonal input -- long winter, but smooth transitions
md.hydrology.englacial_input(1 : end - 1, 🙂 = 3.17e-8 * 3600 * 24 * 365;
for tt = 1 : length(time)
if time(tt) >= 0.4 && time(tt) < 0.7
md.hydrology.englacial_input((1 : end - 1), tt) = -492.75 * cos(2 * pi / 0.3 * (time(tt) - 0.4)) + 492.75 + 1; % m/yr
elseif time(tt) >= 1.4 && time(tt) < 1.7
md.hydrology.englacial_input((1 : end - 1), tt) = -492.75 * cos(2 * pi/0.3 * (time(tt) - 1.4)) + 492.75 + 1; % m/yr
end
end

md.hydrology.moulin_input = zeros(md.mesh.numberofvertices, 1);
%========================================

Based on the constraints on the upper limit of the gap height in the SHAKTI model you described, I guess this must be one of the important reasons why my results are different from yours. As for other reasons, I really don’t know.
I am very excited that SHAKTI can be applied to Antarctica because my research region is in Antarctica. Therefore, I may need your help and answers on model knowledge and operation details in the future, which may delay your time. I hope you can understand! Thanks a lot again!

aleahsommers

Yes, I think the gap height limit may be responsible for the differences. Otherwise, what you have looks pretty good, unless we did something else different in the 2018 example that isn't obvious piecing together the old scripts. Sorry for not having a quick answer to resolve that clearly, but I think you have a good understanding of the setup.

I'll be very happy to help in any way when you'd like to start using SHAKTI in Antarctica. Here is my email address so we can be in touch: aleah.n.sommers@dartmouth.edu

Thanks!