Hi @JohannaBeckmann. I used the PISM friction law when I compared all friction laws in ISSM in this paper:. You should try to initialize with a water layer, for example 2 metres:
md.hydrology.watercolumn_max =2*ones(md.mesh.numberofvertices,1); %m
Feel free to recycle the code I used in my Greenland paper - you will likely need to adjust for your settings.
%Hydrology for PISM friction law
md.hydrology = hydrologypism();
md.hydrology.drainage_rate =1*ones(md.mesh.numberofvertices,1); %mm/a
md.hydrology.watercolumn_max =2*ones(md.mesh.numberofvertices,1); %m
md.initialization.pressure =zeros(md.mesh.numberofvertices,1);
md.initialization.watercolumn =2*ones(md.mesh.numberofvertices,1);
md.basalforcings.groundedice_melting_rate =zeros(md.mesh.numberofvertices,1); %this can be set as a function of geothermal heat or set to some value
%Friction
md.friction=frictionpism();
md.friction.sediment_compressibility_coefficient=0.12*ones(md.mesh.numberofvertices,1);
md.friction.threshold_speed=100; %default = 100 m/a
md.friction.pseudoplasticity_exponent=0.6; %default PISM=0.25; default ISSM after Aschwanden et al 2017: 0.6, 1 means linear sliding law
md.friction.void_ratio=0.69; %default = 0.69 [dimensionless]
md.friction.delta=0.02; %default 0.02 -- lower limit of the effective pressure, expressed as a fraction of overburden pressure [dimensionless]
%Set till friction angle linearly increasing with elevation
%Aschwanden et al 2016 values
% * b_min = -700 m asl
% * b_max = 700 m asl
% * phi_min = 5 degrees
% * phi_max = 40 degrees
% M=(ϕmax−ϕmin)/(bmax−bmin).
% * till_friction_angle = phi_min, where b(x,y) < b_min
% * till_friction_angle = ϕmin+(b(x,y)−bmin) M, where b_min < b(x,y) < b_max (linearly varying depending on bed_elevation)
% * till_friction_angle = phi_max, where b(x,y) > b_max
bed_min = -700; %m asl
bed_max = 700; %m asl
phi_min = 5; %degrees
phi_max = 40; %degrees
M=(phi_max-phi_min)/(bed_max-bed_min); %dependency (slope) of friction angle between bed_min and bed_max
md.friction.till_friction_angle=ones(md.mesh.numberofvertices,1);
%Bed below bed_min m asl
pos=find(md.geometry.bed<bed_min);
md.friction.till_friction_angle(pos)=phi_min;
%Bed above bed_max m asl
pos=find(md.geometry.bed>bed_max);
md.friction.till_friction_angle(pos)=phi_max;
%Bed between bed_min and bed_max
pos=find(md.geometry.bed>bed_min & md.geometry.bed<bed_max);
for pp=1:length(pos)
md.friction.till_friction_angle(pos(pp))=phi_min+(md.geometry.bed(pos(pp))-bed_min)*M;
end
Feel feel to drop me an email if you want more help or would like to discuss the science 🙂 henning.akesson[at]geo.uio.no
Best,
Henning
Reference
Åkesson, H., Morlighem, M., O’Regan, M., & Jakobsson, M. (2021). Future projections of Petermann Glacier under ocean warming depend strongly on friction law. Journal of Geophysical Research: Earth Surface, 126, e2020JF005921. https://doi.org/10.1029/2020JF005921