In my experiment, I'd like to switch the sliding law, from Weertman to Budd, while maintaining velocity and thickness of the glacier. Exponents in the sliding laws are set to 1 (p = q = m = 1) which makes it easy to relate the friction coefficients of the two laws, and the numerical solution does converge within the given number of iterations, meaning the coefficient conversion should be correct. However, after the switch the glacier starts to relax to a new configuration and it shouldn't happen.
It's unclear to me what causes the issue. Any insights? Or better ways to switch sliding laws?
md = loadmodel(org, 'Transient_Steadystate');
md = transientrestart(md);
% current sliding is is Weertman's law
% calculate coefficients in Budd law
C = md.friction.C;
H = md.geometry.thickness;
fric_coef = sqrt(C.^2./(md.materials.rho_ice*md.constants.g.*H));
% enable Budd law
md.friction = friction();
md.friction.p = ones(md.mesh.numberofelements,1);
md.friction.q = ones(md.mesh.numberofelements,1);
md.friction.coefficient = fric_coef;
md.timestepping.final_time = md.timestepping.start_time + 5;
md = solve(md,'tr');