First release of MarshMorpho2D

MarshMorpho2Dexamples.m

@@ -0,0 +1,149 @@
+clear; close all;clc
+%NOTES
+
+%A
+%0: not in the domain
+%1: a normal cell
+%2: the open sea boundary
+
+%Various initiliaztion for plotting%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+cdata = [-1    205  165 0     0     0    0    0          0    255 255 255     1     0    0    0];dlmwrite('mycmap.cpt', cdata, ' ');
+cdata = [-1    205  165 0     0     0    0    0          0    255 255 255     1     0    0    0];dlmwrite('mycmapCLR.cpt', cdata, ' ');
+cdata = [-1    205  165 0     0     0    0    0];dlmwrite('mycmap1.cpt', cdata, ' ');
+
+
+
+
+%%%%%%%%%%%%%%PARAMETERS%%%%%%%%%%%%%%%%%%%%%%%
+P=struct;
+P.kro=0.01; % minimum water depth [m]
+P.DiffS=0.1; %coefficient for tidal dispersion [-]
+
+P.RSLR=1/1000/365;  %Relative Sea Level Rise from mm/yr to m/day (the time unit is the day!)
+
+%tide
+P.Ttide=12.5/24; %[tidal period [day]
+P.Trange=3.1;%1.5; %tidal Trange [m]
+P.TrangeVEG=3.1;%2.7;%1.5; %tidal Trange [m] MOST OF THE TIME THIS IS THE SAME OF TIDAL RANGE
+
+%SSC at the sea boundary
+P.co2=5/1000; %Sea boundary SSC for mud [g/l]
+
+P.ws2=0.2/1000;%mud settlign velocity [m/s]
+P.rbulk2=800; %dry bulk density [kg/m3]
+
+%mud parameters
+P.me=10^-5*24*3600;  %mud erodability kg/m2/day!!!
+P.tcrgradeint=0.2; %linear increase in tcr below MLW [pa/m]
+P.taucr=0.2; %crtical shear stress unvegetated[Pa]
+P.crMUD=3.65/365; %[m2/day] unvegetated creep coefficient
+P.crMARSH=0.2/365; %[m2/day] vegetated creep coefficient
+
+%Vegetation parameters
+P.dBlo=-(0.237*P.TrangeVEG-0.092);%lower limit for veg growth [m]  see McKee, K.L., Patrick, W.H., Jr., 1988. 
+P.dBup=P.TrangeVEG/2;%upper limit for veg growth [m] ee McKee, K.L., Patrick, W.H., Jr., 1988. 
+P.Cb=0.02;% manning coefficient unvegegated
+P.Cv=0.1; %mannign coefficenitn in vegegated areas
+P.wsB=1/1000;%settlign velocity in vegetated area [m/s]
+P.taucrVEG=0.5;%Critical stress for vegetated areas
+
+P.Korg=8/1000/365; %organci sediment production [m/day]
+
+%to impose the boundary condition. 
+%if =1 the boundary cell becomes equal to the closest cell in the domain 
+%if = 0 the boundary depth does not change
+P.imposeseaboundarydepthmorphoALL=1;
+
+tmax=500;%how many time steps in total
+dtO=2*365;%length of one time step
+tINT=1;% how many time step you shoudl plot
+time=[1:tmax]*dtO/365;
+
+%to create a video (no or yes)
+makevideo=0;
+
+
+%%%%%%%%%%%%%%%Geometry Initilization%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%here are two examples. chose one or make your own grid
+[N,M,dx,A,z,x,y,msl]=initializegeometrySWENEY(P);
+%[N,M,dx,A,z,x,y,msl]=initializegeometryWEST(P);
+
+
+IO.z=z;
+IO.msl=msl;
+fIO.FLX=0;%keep track of sediment flux through open boundary
+fIO.KBTOT=0;%keep track of how much organic was produced
+
+
+zoriginal=z;%store the original bed elevation
+sumY2IN=sum(z(A==1));%Store value for mass balance check
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+%%%MAIN LOOP
+if makevideo==1;v=VideoWriter('PUTHEREVIDEONAME','Motion JPEG AVI');open(v);end
+s=0;
+for t=1:tmax;  
+
+%THIS IS THE CORE OF THE MODEL!!!!!!!!!
+[IO,fIO,PLT]=mainevolutionstep(A,P,dx,dtO,IO,fIO);
+
+
+
+
+
+
+%%%PLOT
+if mod(t,tINT)==0;s=s+1;  
+
+%read the variables
+names = fieldnames(IO);
+for i=1:length(names);eval([names{i} '=IO.' names{i} ';' ]);end
+
+%read the fluxes
+names = fieldnames(fIO);
+for i=1:length(names);eval([names{i} '=fIO.' names{i} ';' ]);end
+
+%read the plot
+names = fieldnames(PLT);
+for i=1:length(names);eval([names{i} '=PLT.' names{i} ';' ]);end
+
+
+ax1 = subplot(1,3,1); 
+IM=imagesc(y,x,-(z+msl));axis equal;set(IM,'alphadata',~(A==0));set(gca,'YDir','normal');%colormap('jet');
+cmp=demcmap([-3 P.Trange/2],256); %-3 1
+colormap(ax1,cmp)
+caxis([-3 P.Trange/2]);
+title(strcat(num2str(time(t)),' years'))
+
+ax1 = subplot(1,3,2);
+IM=imagesc(y,x,-(zoriginal));axis equal;set(IM,'alphadata',~(A==0));set(gca,'YDir','normal');%colormap('jet');
+cmp=demcmap([-3 P.Trange/2],256);
+colormap(ax1,cmp)
+caxis([-3 P.Trange/2]);
+title('original bathymetry')
+
+ax2 = subplot(1,3,3);
+IM=imagesc(y,x,1000*SSC);set(IM,'alphadata',~(A==0));axis equal;set(gca,'YDir','normal');%colormap('jet');%colorbar('hori') 
+caxis([0 20]);colormap('jet')
+colormap(ax2,flipud(hot))
+title('SSC [mg/l]')
+
+
+sumY2=sum(z(A==1));
+sumFLUX2=-FLX*dx;
+%NOTE: Thsi is the equivalent volumetric flux, not the mass flux
+checksum=[-(sumY2IN-sumY2)+sumFLUX2/dx^2]+KBTOT
+%if close to zero it means there are no issues (mass is conserved!!!)
+
+if makevideo==1;V=getframe(figure(1));writeVideo(v,V);end
+pause(0.1)
+end
+
+end
+
+if makevideo==1;close(v);end %UNCOMMENT THIS TO CREATE A VIDEO
+
+
+

bedcreep.m

@@ -0,0 +1,36 @@
+function z=bedcreep(z,A,crMUD,crMARSH,dx,dt,VEG);
+
+
+creep=A*0;
+creep(VEG==0)=crMUD;           
+creep(VEG==1)=crMARSH;
+
+D=(creep)/(dx^2)*dt;
+
+G=0*z;
+p=find(A==1);%exclude the NOLAND CELLS (A==0)
+NN=length(p);G(p)=[1:NN];rhs=z(p);[N,M]=size(G);i=[];j=[];s=[];
+
+S=0*A;
+[row col]=ind2sub(size(A),p);
+for k = [N -1 1 -N] 
+
+[a,q]=excludeboundarycell(k,N,M,p);
+a=a(A(q(a))==1);%only inclued the cells in whcih you can creep to
+
+
+gradF=abs(z(p(a))-z(q(a)))/dx;
+facNL=gradF>=tan(1/180*pi);
+
+value=(D(p(a))+D(q(a)))/2.*facNL;
+
+S(p(a))=S(p(a))+value; %exit from that cell
+i=[i;G(q(a))]; j=[j;G(p(a))]; s=[s;-value]; %gain from the neigborh cell
+end
+
+%summary of the material that exits the cell
+i=[i;G(p)]; j=[j;G(p)]; s=[s;1+S(p)];
+ds2 = sparse(i,j,s);%solve the matrix inversion
+P=ds2\rhs;z(G>0)=full(P(G(G>0)));
+
+

excludeboundarycell.m

@@ -0,0 +1,8 @@
+function [a,q]=excludeboundarycell(k,N,M,p);
+[row col]=ind2sub([N M],p);
+if k==N;a=find(col+1<=M);end;
+if k==-N;a=find(col-1>0);end;
+if k==-1;a=find(row-1>0);end;
+if k==1;a=find(row+1<=N);end;
+
+q=p+k;%the translated cell

getwaterdepth.m

@@ -0,0 +1,18 @@
+function [h,ho,fTide,dtide,dHW,wl]=getwaterdepth(range,msl,z,kro);
+DH=range; %total water displacement
+dHW=max(0,z+msl+range/2);%water depth at MHW
+dtide=max(0,z+msl+range/2);%water depth at MHW
+fTide=min(1,max(10^-3,dHW/DH));%hydroperiod
+h=0.5*(max(0,dHW)+max(0,dHW-DH)); %for the flow
+%ho=h; %the original water depth without the limiter
+h(h<kro)=kro;%reference water depth for flow calculation
+ho=h;
+
+%the depth that gives the tidal prism
+ho(-z>(msl+range/2))=0;
+
+wl=h-z;
+
+% figure;
+% imagesc(ho)
+% pause

initializegeometrySWENEY.m

@@ -0,0 +1,19 @@
+function [N,M,dx,A,z,x,y,msl]=initializegeometry(P);
+
+dx=2;
+load SWENEY;
+z=-zSW3(1:end-10,:);z=double(z);
+
+%attach
+z(end-25:end,132:141)=-1.35;
+
+[N,M]=size(z);
+A=ones(N,M);
+A(z==-999 | z<-1.9)=0;
+A(end,141:155)=2;
+A(isnan(z))=0;
+z(A==0)=-2;
+
+x=[0:N-1]*dx/1000;y=[0:M-1]*dx/1000;
+
+msl=0;

initializegeometryWEST.m

@@ -0,0 +1,27 @@
+function [N,M,dx,A,z,x,y,msl]=initializegeometry(P);
+
+load WEST;
+z=-PIE2mBELLO;
+z(247:249,361:363)=-1.6;%tah house
+z=z(:,1:end-24);
+
+
+dx=2;
+[N,M]=size(z);
+
+%%%%%%%%%%%cell types
+A=ones(N,M);
+
+A(z==-999)=0;
+A(:,end)=0;
+A(95*4/dx:114*4/dx,end)=2;
+
+A(isnan(z))=0;
+
+z(A==0)=-2;
+
+
+x=[0:N-1]*dx/1000;y=[0:M-1]*dx/1000;
+
+
+msl=0;

mycmap.cpt

@@ -0,0 +1 @@
+-1 205 165 0 0 0 0 0 0 255 255 255 1 0 0 0

mycmap1.cpt

@@ -0,0 +1 @@
+-1 205 165 0 0 0 0 0

mycmapCLR.cpt

@@ -0,0 +1 @@
+-1 205 165 0 0 0 0 0 0 255 255 255 1 0 0 0

seaboundaryNeumanbedelevationALLBOUNDARY.m

@@ -0,0 +1,15 @@
+function [z]=seaboundaryNeumanbedelevationALLBOUNDARY(A,z)
+
+[N,M]=size(A);
+p=find(A==2);%exclude the NOLAND CELLS (A==0)
+
+[row col]=ind2sub(size(A),p);
+for k = [N -1 1 -N] 
+
+[a,q]=excludeboundarycell(k,N,M,p);
+a=a(A(q(a))==1);%only inclued the cells in whcih you can creep to
+
+z(p(a))=z(q(a));
+end
+
+

sedtran.m

@@ -0,0 +1,94 @@
+function [EmD,SSM,FLX]=sedtran(d,A,DiffS,h,ho,E,WS,dx,dt,rbulk,co,Ux,Uy,FLX,fTide,Ttide,kro);
+
+A(ho<=0)=0; %eliminate the cells in which the water depth is too small
+
+%rivermouthfront
+p = find(A>0);%exclude the NOLAND CELLS (A==0)
+G=0*d;NN=length(p);G(p)=[1:NN];
+rhs=E(p); %in the rhs there is already the additon of the erosion input
+[N,M]=size(G);
+i=[];j=[];s=[];S=0*G;
+
+%boundary conditions imposed SSC
+a=find(A==2);rhs(G(a))=co.*h(a).*fTide(a);
+
+Dxx=(DiffS*Ttide/2*(abs(Ux.*Ux))*(24*3600).^2)/(dx^2).*h;%.*(ho>kro);%.*(hgross>0.01);%% the h is not the coefficient for diffusion, it the h in the mass balance eq.
+Dyy=(DiffS*Ttide/2*(abs(Uy.*Uy))*(24*3600).^2)/(dx^2).*h;%.*(ho>kro);%.*(hgross>0.01);
+
+%the factor 24*3600 is used to convert the Ux and Uy from m/s to m/day
+
+[row col]=ind2sub(size(A),p);
+for k = [N -1 1 -N]  %go over the 4 directions for the gradients
+
+%avoid to the the cells out of the domain (risk to make it periodic...)
+[a,q]=excludeboundarycell(k,N,M,p);
+
+a=a(A(q(a))==1 | A(q(a))==2 | A(q(a))==10);%exlcude the translated cell that are NOLAND cells
+
+%go over the extradiagonal direction for each gradient direction
+if (k==N | k==-N);D=Dyy;else;D=Dxx;end;
+
+DD=(D(p(a))+D(q(a)))/2;
+
+Fin=0*DD;Fin(A(p(a))==1)=1; % (A(q(a))==1 | A(q(a))==10)=1; to conserve the mass a the river mouth = no input
+Fout=0*DD;Fout(A(q(a))==1)=1; %needed not to affect the b.c. -> Do not choose 2 and 1p
+
+%tidal dispersion component
+value=DD./h(p(a))./fTide(p(a));
+
+S(p(a))=S(p(a))+value.*Fin; %exit from that cell
+i=[i;G(q(a))]; j=[j;G(p(a))]; s=[s;-value.*Fout]; %gain from the neigborh cell
+
+end
+
+%summary of the material that exits the cell
+settling=24*3600*WS(p)./h(p).*fTide(p);%.*(h(p)>3);
+
+%sea boundary
+a=find(A(p)==2);%find the co b.c.
+settling(a)=0;%do not settle in the b.c. (to not change the SSC)
+S(p(a))=1;%to impose the b.c.
+
+
+i=[i;G(p)]; j=[j;G(p)]; s=[s;S(p)+settling];
+ds2 = sparse(i,j,s);%solve the matrix inversion
+P=ds2\rhs;%P=mldivide(ds2,rhs);
+
+SSM=0*d;SSM(G>0)=full(P(G(G>0)));%rescale the matrix
+
+%update the bed
+EmD=0*A;EmD(p)=(E(p)-SSM(p).*settling)/rbulk;
+
+
+
+
+
+
+
+
+
+%%%%%%%%OUTPUT FOR SSM BALANCE. DOES NOT AFFECT COMPUTATION!!!
+p = find(A==1);[row col]=ind2sub(size(A),p);
+
+%sea boundary
+Q=0;
+for k = [N -1 1 -N]
+
+[a,q]=excludeboundarycell(k,N,M,p);
+
+%only get the cell that discharge into a b.c. A==2
+a=a(A(q(a))==2);
+
+%for each gradient direction
+if (k==N | k==-N);D=Dyy;else;D=Dxx;end;
+DD=(D(p(a))+D(q(a)))/2;
+
+Q=Q+sum(DD.*(SSM(p(a))./h(p(a))./fTide(p(a))-SSM(q(a))./h(q(a))./fTide(q(a)))); %exit from that cell
+
+end
+
+FLX=FLX+dx*dt*Q/rbulk; %(note the the divided dx is for the gradient, not for the cell width!)
+%%%%%%%%%%%%%%%%%%%
+
+
+