%--------------------------------------------------------------------------------------
%Large–scale volcanic deposit fluidisation by dilute pyroclastic density currents
%by Karim Kelfoun & Antonio Proaño
%Nature Geoscience
%Done by K.Kelfoun, LMV-UCA, 2023
%--------------------------------------------------------------------------------------

%This code is made for the simulation of gas diffusion
%using an upwind scheme of finite differences
%For clarity, the code has been commented and simplified
%to keep only the essentials of the original code

%Runs in Matlab https://fr.mathworks.com/products/matlab.html 
%     or Octave https://octave.org/

clear all
close all

%Granular bed characteristics
H  = 0.2;               %thickness of the granular bed (m)
dz = 0.001;             %mesh size (m)
z  = dz/2:dz:H+dz;      %elevation of each mesh
rhog = 1450;            %density of the granular bed (kg/m3)
epsilon = z*0+0.3;      %connected porosity
epsilon(end)=1;         %the last porosity is 100% = atmosphere over the bed

mu = 1.85e-5;           %air viscosity
K  = 1e-12;             %gas permeability in the granular bed

%Parameters of the surge decompression
Dt   = 10;              %time of surge decompression (s)
Patm = 101300;          %atmospheric pressure (Pa)
Pini = Patm + 20000;    %initial pressure (Ps0 in the paper) of the base of the surge (Pa)

%Initial conditions
P = z*0+Pini;           %initial gas pressure for all meshes
rho = P/Patm*1.2;       %air density related to the pressure (from the ideal gas law)

g=9.81;                 %gravity

%Time parameters of the simulation
tmax = 15;             %duration of the simulation
dt = 1e-5;             %time step - to be adjusted if you change the bed geometry

M(1:length(z), 1:round(tmax*100)) = 0; %Matrix used to record the pressure calculated

%--------- Simulation of the depressurisation with time -------------------
for t = dt:dt:tmax

    %Mean pressure, mean porosity and mean density at the cell edges
    %where are calculated fluxes and velocities
    Pm = ( P(2:end)+P(1:end-1) )/2;
    epsilonm = ( epsilon(2:end)+epsilon(1:end-1) )/2;
    epsilonm(1) = epsilonm(2);
    rhom = ( rho(2:end)+rho(1:end-1) ) / 2;

    %Velocity of the gas in the pores
    u = K./epsilonm   ./mu.*(  ( P(1:end-1)-P(2:end) )./dz  -  rhom*g    );  %Eq. 2 of the paper

    %Upstream scheme for the fluxes
    drho_dt = Pm*0;
    ind = find(u>0); %if the gas moves upwards
    drho_dt(ind) = u(ind).*rho(ind).*epsilonm(ind) ;   %mass fluxes, in kg/s/m²
    ind = find(u<0); %if the gas moves downwards
    drho_dt(ind) = u(ind).*rho(ind+1).*epsilonm(ind) ; %mass fluxes, in kg/s/m²

    %New densities in the cells (Eq. 1)
    rho(2:end-1) = rho(2:end-1) + (drho_dt(1:end-1)-drho_dt(2:end))./dz*dt./epsilon(2:end-1); %kg/m3
    rho(1)       = rho(1)       - drho_dt(1)   ./dz*dt./epsilon(1);                           %kg/m3
    rho(end)     = rho(end)     + drho_dt(end) ./dz*dt./epsilon(end);                         %kg/m3

    %New pressures in the cells (Eq. 3)
    P   = rho.*Patm/1.2;

    %Boundary condition: pressure over the granular bed
    P(end) = (Pini-Patm)*( 1-t/Dt )+Patm; %for the pressure (Ps + Patm in the paper)
    if P(end)<Patm
        P(end)=Patm;
    end
    rho(end)=P(end)/Patm*1.2; %... and the density
    %----------------------------------------------------------------------

    if abs(t*1-round(t*1))<1e-6 %time step for plotting
        plot(P, z)
        hold on
        fprintf('%3.0f\n', t)
        xlabel('Pressure (Pa)')
        ylabel('Elevation (m)')
        drawnow
    end
    if abs(t*100-round(t*100))<1e-6 %time step for recording
        M(:, round(t*100)) = P';
    end
end
%--------------- End of the simulation loop -------------------------------

%Figure 1 - Pressure
figure(1)
imagesc(100*dt:100*dt:tmax, z, M);
set(gca, 'ydir', 'normal')
xlabel('time (s)')
ylabel('elevation (m)')
title('Pressure in the granular bed (Pa)')
colorbar

%Figure 2 - Friction angle
figure(2)
h=H-z;                                          %depth
Psolid = rhog*9.81*h;                           %solid pressure
Psolid = Psolid'*(1+[1:size(M,2)]*0);           %solid pressure - size of M
dM=M-(1+[1:size(M,1)]*0)'*M(end,:);             %Pressure difference P-Ps
phi=atand(tand(27)*(Psolid-dM)./Psolid);        %friction angle
phi(phi<0)=0;
imagesc(1000*dt:1000*dt:tmax, z, phi)
set(gca, 'ydir', 'normal')
title('Friction angle of the granular bed, °')
xlabel('time (s)')
ylabel('elevation (m)')
colorbar
hold on
contour(1000*dt:1000*dt:tmax, z, phi, [1e-6 1e-6], 'w');

