# Unverified solution in MATLAB

From Interactive System for Ice sheet Simulation

% Possible solution to Kees' problem in MATLAB % % Based on eqns. presented on % http://websrv.cs.umt.edu/isis/index.php/Kees'_assignment % Possibly useful as a solution template for those familiar with MATLAB. % Note possible problem in lines 55-58. clear all close all clc % Define constants. g = 9.8; % m/ s^ 2; gravitational acceleration rho = 920; % kg/ m^ 3; density of ice delta_x = 100; % m; distance step x_max = 25* 10^ 3; % m; length of problem domain delta_t = 10^ -2; % yr; time step (needs to be small, about 10^ -2 yr) t_max = 10^ 1; % yr; integration time of model M_0 = 4; % m; net accumulation at upper end of problem domain S_M = 4/ 10^ 4; % 1/ m; net accumulation gradient (note change from % value given on wiki page, above) bed_slope = 0.1; % slope of glacier bed Glen_n = 3; % d'less; exponent in Glen's flow law Glen_A = 10^ -16; % Pa^ -3/ yr; coefficient in Glen's flow law % Calculate derived constants. C = (2/ (Glen_n+ 2))* Glen_A* (rho* g)^ Glen_n; % Set up initial condition and plotting variables. Initialize vectors. x = 0: delta_x: x_max; % m; distances from upper end of domain for % each node t = 0: delta_t: t_max; % yr; times corresponding to each time step z_bed = wrev(bed_slope* x); % m; bed height as a fn of x num_x = numel(x); % number of nodes num_t = numel(t); % number of time steps D = zeros(1, num_x- 1); % m^ 2/ yr; diffusivities at staggered points H = zeros(1, num_x); % m; thickness of ice at different points M = zeros(1, num_x); % m/ time step; net mass balance per time step num_D = numel(D); % number of points in offset grid phi = zeros(1, num_D); % m^ 2/ yr?; flux of matter between nodes % Calculate the net mass balance per time step within the problem domain. M = (M_0- S_M* x)* delta_t; % m/ time step % Step through time. for count1 = 1: 1: num_t; for count2 = 1: 1: num_D; D(count2) = C* ((H(count2)+ H(count2+ 1))/ 2)^ (Glen_n+ 2)* ... ((H(count2+ 1)+ z_bed(count2+ 1))- ... (H(count2)+ z_bed(count2)))^ (Glen_n- 1); phi(count2) = D(count2)* (((H(count2+ 1)+ z_bed(count2+ 1))- ... (H(count2)+ z_bed(count2)))/ delta_x); end for count3 = 2: 1: (num_x- 1); H(count3) = H(count3)+ (phi(count3)- phi(count3- 1))/ ... delta_x+ M(count3); % H(count3) = H(count3)- (phi(count3)- phi(count3- 1))/ ... % delta_x+ M(count3); end H(1) = 0; for count4 = 1: 1: num_x; if H(count4) < 0; H(count4) = 0; end end end % Plot results. ice_top = z_bed+ H; % m; elevation of ice surface above lowest node plot(x, ice_top, 'b-') hold on plot(x, z_bed, 'k-') xlabel('Horizontal distance (m)') ylabel('Elevation (m)')