%Esther Adelstein, EPS 109 Fall 2012 Final Project
%Modeling Post-Glacial Rebound for Varying Depression Wavelengths
%In order to use linear elastic relationship, assume W<<lambda
%vertical depression = positive W, negative pressure. Return to
%undeformed W = 0 approximated by viscous flow in the half-space
%Assume mantle is an incompressible (const. density)/Newetonian fluid, so use continuity
%equation: through an infinitesimal rectangular element,
%d(vx)/dx+d(vy)/dy=0. 
%Stream function satisfies the above, takes form of
%sin(2*pi*x/lambda)*Y(y), plug into biharmonic eqn, (del)^4(streamfunction)
%= 0...integrate to solve for stream function, differentiate to find
%velocity components. Important here is equating hydrostatic pressure head
%due to topography to normal stress at upper boundary of fluid half-space,
%-rho*g*W=p-2*mu*dvdy. Here and for some boundary conditions, W<<lambda
%assumption is necessary
%Assume intial displacement W = W0*cos(2*pi*x/lambda)
%W(t) = W*exp(-lambda*rho*g*t/(4*pi*mu))
%Reference: Turcotte,Schubert, Geodynamics second edition


clear
clf

% define variables
mu = 1e21; % Pa*s, typical mantle viscosity
rho = 3300; % kg/m^3, typical mantle rock density
g = 9.8;  % m/s^2, gravitiational acceleration
lambda = 5e6; % m, estimated wavelength of depression from Laurentide ice sheet

% set initial conditions
N = 100;
x = [0:lambda/N:lambda];
y = [0:lambda/N:lambda];
W = zeros(N+1,N+1);
W0 = 61; % m, initial vertical depression 
for j = 1:N+1
    for i = 1:N+1
        W(i,j) = (W0*cos(2*pi*i/(N+1))-W0)/2;% modified from expression derived above to set W = 0 at surface
        if j<=(i-(N+1)/2)^2/((N+1)/4)% parabolic boundary of depression
            W(i,j) = 0;
        end
        lambdai(i,j) = abs(2*sqrt(j)/(N+1)*lambda);% set wavelength for each cell equal to the width of the parabola at the row where the cell resides;
    end
end

% iterate over time
f = 0; % counter for time
g = 0; % counter for images
dt = 100; %years
for t = 0:dt:1e4; % assume Laurentide ice sheet melted 10,000 years ago
    for i = 1:N+1
        for j = 1:N+1
            W(i,j) = W(i,j)*exp(-lambdai(i,j)*rho*g*dt*1e7/(4*pi*mu));% W updated after each timestep, so use dt here
        end
    end
    % surface plot
    xx = 0:lambda/N:lambda;
    yy = 0:lambda/N:lambda;
    [X,Y] = meshgrid(yy,xx);
    subplot(1,2,1);
    surf(X,Y,W);
    view([+50 10]);
    xlabel('horizontal location (m)')
    ylabel('horizontal location (m)')
    zlabel('vertical depression (m)')
    title(sprintf('Post-Glacial Rebound, t = %5.0f years',t));
    axis([0 lambda 0 lambda -80 0])
    pause(0.1)
    % topography plot
    subplot(1,2,2);
    pcolor(X,Y,W)
    shading interp
    contour(X,Y,W,-200:10:100,'k')
    [c,h] = contour(X,Y,W,-200:10:100,'k');
    clabel(c,h)
    xlabel('horizontal location (m)')
    ylabel('horizontal location (m)')
    zlabel('vertical displacement (m)')
    title('Topography')
    pause(0.01)
    % generate images
    if t == 0 
       saveas(gcf,['post_glacial_rebound',sprintf('%03d',g),'.png'])
    end
    if t == 1e4 
       saveas(gcf,['post_glacial_rebound',sprintf('%03d',g),'.png'])
    end
    if mod(t+5e2,1e3) == 0
       saveas(gcf,['post_glacial_rebound',sprintf('%03d',g),'.png'])
    end
    g = g+1;
    f = f+dt;
end
% This violates the W<<lambda assumption  for smaller lambda, but results are reasonable for larger lambda  