# Coastline evolution #------------------------------------------------------------------------------ rhos=2650 # Sediment mass density = 2650 kg/m3 rho=1025 # Sea water mass density = 1025 kg/m3 g=9.806 # Gravity acceleration = 9.806 m/(s2) n=0.4 # Porosity n=0.6 K=0.7 # Dimensionless parameter to compute the actual sediment transport s=1/30 #Beach slope; it is not used, unless one uses other formulas than CERC d50m=0.0005 #Beach grain diameter; it is not used, unless one uses other formulas than CERC # Time conversion day/s = 86400 t=0 pr1=scan("coast.txt") #Upload of the coastline geometry pr1=array(pr1,dim=c(2,20)) dy=scan("wave-forcing.txt",n=1) #Uplod of dy dt=scan("wave-forcing.txt",skip=1,n=1) #Upload of dt ndeltat=scan("wave-forcing.txt",skip=2,n=1) #Upload of the number of time steps of the simulation pr2=scan("wave-forcing.txt",skip=3) #Upload of wave forcing pr2=array(pr2,dim=c(5,20)) xlc=pr1[1,] #position of coastline xlcl=pr1[2,] #position of limit coastline #xlc=cumsum(rep(dy,20)) #Test instruction - allows users to upload a straight coastline with a slope of 45 degrees dlf=pr2[1,] #Position of breakerline Hb=pr2[2,] #wave height at the breakerline lf=pr2[3,] #wave length at the breakerline db=pr2[4,] #sea bottom depth at the breakerline alpha=pr2[5,] #wave angle with respect to y - positive when waves are oriented towards increasing Ys y=cumsum(rep(dy,20)) xlcold=xlc #Instruction to keep into memory the original position of the shoreline plot(xlcold,y,type="l",ylim=c(0,3000),xlim=c(350,1500)) #Plot instructions - users can remove them lines(xlcl,y,col="red") ql=rep(0,20) for (j in 1:19) { xlc[j] = (xlc[j]+xlc[j+1])/2 } for(i in 1:ndeltat) #Cycle along the timesteps - set the number of iteration - suggested limits are 1:ndeltat { for (ii in 2:19) #Cycle along the spatial steps - along the shoreline { xlcaux=xlc[ii+1]-xlc[ii-1] alpha2=alpha[ii]-atan(xlcaux/2/dy) #Computation of the wave angle with respect to average coastline Hbaux=(Hb[ii-1]+Hb[ii]+Hb[ii+1])/3 #Computation of average wave height #alpha2=pi/2-pi/4 #Test instruction alpha=rep(pi/2+p/4,20) - allows users to set a specific wave angle with respect to shoreline - setting alpha2=pi/2 should lead to no changes in the shoreline morphology Eb=1/8*rho*g*Hbaux^2 Ebw=Eb*(g*db[ii])^0.5 #ql[ii]= (1.28*s/d50m)*(Hbaux**(3.5))*sin(2*alpha2)/86400*dt #Computation of volumetric sediment transport with Kamphuis formula ql[ii]=(K*0.5*Ebw*sin(2*alpha2))/((rhos-rho)*g*(1-n)) #Computation of volumetric sediment transport rate with CERC formula } #ql is positive if alpha2 is lower than pi/2. In this case the sediment transport is from right to left in the given exercise. It is negative if alpha2 is greater than pi/2, and then the transport is from left to right dql=diff(ql) #Be careful, the resulting vector is counting 1 less element with respect to ql #cat(ql,fill=T) #Test instruction to print the volumetric sediment tranport rate #cat(dql,fill=T) #Test instruction to pring the differenced volumetric sediment transport rate for (j in 1:19) { dcsi=-1/db[j]*dt*dql[j]/dy #The negative sign is to ensure coherence. If alpha2 is less than pi/2, the sediment rate is positive. Therefore, a negative increment means that there is deposition and therefore accretion of the beach. #Conversely, a positive increment means that there is erosion and therefore erosion of the beach. #The area of the parallelogram is computed here xlc[j] = xlc[j]+dcsi } #enddo } #enddo lines(lowess(xlc[1:19],y[1:19]+dy/2,f=0.38),col="blue")