


GETALPHA2 - compute the friction parameter alpha2
This routine calculates the basal friction coefficient
o alpha2= drag^2 * Neff ^r * vel ^s
with
o Neff=rho_ice*g*thickness+rho_ice*g*bed,
o r=q/p and s=1/p
Usage:
alpha2=Getalpha2(parameters)

0001 function alpha2=Getalpha2(parameters) 0002 %GETALPHA2 - compute the friction parameter alpha2 0003 % 0004 % This routine calculates the basal friction coefficient 0005 % o alpha2= drag^2 * Neff ^r * vel ^s 0006 % with 0007 % o Neff=rho_ice*g*thickness+rho_ice*g*bed, 0008 % o r=q/p and s=1/p 0009 % 0010 % Usage: 0011 % alpha2=Getalpha2(parameters) 0012 0013 %recover parameters 0014 gravity=parameters.g; 0015 rho_ice=parameters.rho_ice; 0016 rho_water=parameters.rho_water; 0017 K=parameters.k; 0018 bed=parameters.b; 0019 thickness=parameters.h; 0020 velocity_x=parameters.velocities(:,1); 0021 velocity_y=parameters.velocities(:,2); 0022 if size(parameters.velocities,2)==3, %we have to consider vz if Stokes model 0023 velocity_z=parameters.velocities(:,3); 0024 end 0025 0026 %From bed and thickness, compute effective pressure when drag is viscous: 0027 Neff=gravity*(rho_ice*thickness+rho_water*bed); 0028 0029 %If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 0030 %the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 0031 %for this should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 0032 %replace it by Neff=0 (ie, equival it to an ice shelf)*/ 0033 pos=find(Neff<0); 0034 Neff(pos)=0; 0035 0036 %recover p and q in the drag law u = k * sigma ^ p * Neff ^ q : 0037 pcoeff=parameters.p; 0038 qcoeff=parameters.q; 0039 rcoeff=qcoeff./pcoeff; 0040 scoeff=1./pcoeff; 0041 0042 %We need the velocity magnitude to evaluate the basal stress: 0043 if size(parameters.velocities,2)==2, 0044 velocity_mag=sqrt(velocity_x.^2+velocity_y.^2); 0045 elseif size(parameters.velocities,2)==3, 0046 velocity_mag=sqrt(velocity_x.^2+velocity_y.^2+velocity_z.^2); 0047 end 0048 0049 alpha2=(K.^2).*(Neff.^rcoeff).*velocity_mag.^(scoeff-1);