function [QS]=transport(S,dPHI,HS,varargin) % function [QS]=transport(S,dPHI,HS,xp,yp) % % INPUT : % S : Structure with input information of the ShorelineS model, with relevant fields: % .dnearshore: depth at toe of dynamic profile % .trform : transport formulation (either 'CERC', 'KAMP', 'MILH', 'CERC3', 'VR14') % .b : Calbration factor of CERC formula (only CERC) % .tper : Wave period [s] % .d50 : Median grain size [um] % .tanbeta : Slope [1:n] (only Kamphuis and Milhomens) % .Pswell : Relative part of the Pswell (in percentage, specify a value between 0 and 100) % .rhos : Density of sediment [kg/m3] % .rhow : Density of water [kg/m3] % .gamma : Breaking coefficient [-] % .porosity : Porosity [-] % .g : Acceleration of gravity [m/s2] % dPHI : Relative angle of offshore waves with respect to the coastline (corrected to 0 in shadow zones) % HS : Breaking wave height (or diffracted wave height) [Nx1] % dPHI : Relative angle of breaking waves (or diffracted wave angle) % OUTPUT : % QS : Transport rates in grid cells [1xN] (in [m3/yr] including pores) % % Copyright : IHE-delft & Deltares, 2020-feb % Authors : Bas Huisman, Dano Roelvink % License : GNU GPLv2.1 %% Compute gradient of HS % why can this not be computed elsewhere? % e.g. in the angles function? % Why do we need a seperate HS and dPHI? if S.wave_diffraction==1 && ~isempty(varargin) xp = varargin{1}; yp = varargin{2}; n=length(HS); for i=1:length(HS) im1=max(i-1,1); ip1=min(i+1,n); dHS(i)=(HS(ip1)-HS(im1))/( hypot(yp(ip1)-yp(im1),xp(ip1)-xp(im1)) ); end end %% Transport : CERC with the offshore wave height and direction if strcmpi(S.trform,'CERC') k=0.2; % using CERC (1984) value of k(SPM,Hs)=0.39 is suggested, but this is typically quite high QS = S.qscal .* S.b .* HS.^2.5 .* sin(2*dPHI); % this formulation uses the offshore wave height and direction if S.wave_diffraction==1 % use HS and dPHI, and second-order component dHS QS = S.qscal .* S.b .* HS.^2.5 .* (sin(2*dPHI)-2*cos(dPHI).* dHS); end %% Transport : CERC with the offshore wave height and direction, including an implicit refraction from offshore to nearshore within the transport formulation elseif strcmpi(S.trform,'CERC2') k=0.39; % using CERC (1984) value of k(SPM,Hs)=0.39 b1 = k .* (S.rhow .* S.g.^0.5 ./ (16 .* sqrt(S.gamma).* (S.rhos-S.rhow) .* (1-S.porosity))); % = theoretical 'b1' factor from Shore Protection Manual b2 = b1 .* ((S.gamma.*S.g).^0.5 ./(2*pi)).^0.2; % b_theory = 0.0946 * 365*24*60*60 = 2.9833E+6 QS = S.qscal .* 365*24*60*60*b2.*HS.^(12/5).*S.tper.^(1/5).*cos(dPHI).^(6/5).*sin(dPHI); %% Transport: CERC with the nearshore breaking wave height and direction (e.g. Vitousek, Barnard, 2015) elseif strcmpi(S.trform,'CERC3') QS = S.qscal .* 365*24*60*60 .* 0.023 .* S.g.^0.5 .* (S.gamma).^-0.52 .* HS.^2.5 .* sin(2.*dPHI); %% Transport : Kamphuis elseif strcmpi(S.trform,'KAMP') QSkampmass=S.qscal .* 2.33 .* S.tper.^1.5 .* S.tanbeta.^0.75 .* S.d50.^-0.25 .* HS.^2 .* (abs(sin(2.*dPHI)).^0.6.*sign(dPHI)); % <- this is the immersed mass under water QS = 365*24*60*60*(QSkampmass ./(S.rhos-S.rhow)) ./(1.0-S.porosity); if S.wave_diffraction==1 % use HS and dPHI, and second-order component QSkampmass=S.qscal .* 2.33 .* S.tper.^1.5 .* S.tanbeta.^0.75 .* S.d50.^-0.25 .* HS.^2 .* ( abs(sin(2*dPHI)).^0.6.*sign(dPHI) - (2/S.tanbeta).*cos(dPHI).*dHS ); QS = S.qscal .* 365*24*60*60*(QSkampmass(i) /(S.rhos-S.rhow)) /(1.0-S.porosity); end %% Transport : Mil-Homens (2013) elseif strcmpi(S.trform,'MILH') QSmilhmass=S.qscal .* 0.15.* S.tper.^0.89 .* S.tanbeta.^0.86 .* S.d50.^-0.69 .* HS.^2.75 .* (abs(sin(2.*dPHI)).^0.5.*sign(dPHI)); % <- this is the immersed mass under water QS = 365*24*60*60*(QSmilhmass ./(S.rhos-S.rhow)) ./(1.0-S.porosity); %% Transport : Van Rijn (2014) elseif strcmpi(S.trform,'VR14') vtide=0; kswell=0.015.*S.Pswell+(1-0.01.*S.Pswell); vwave=0.3.*real(sin(2.*dPHI)).*(S.g.*HS).^0.5; vtotal=vwave+vtide; QSvr14mass=S.qscal .* 0.0006 .* kswell .* S.rhos .* S.tanbeta.^0.4 .*S.d50.^-0.6 .* HS.^2.6 .* abs(vtotal).*sign(dPHI); QS = 365*24*60*60*(QSvr14mass ./(S.rhos-S.rhow)) ./(1.0-S.porosity); end %% Set transport to 0 when angle exceeds 180 degrees QS(abs(dPHI)>pi/2)=0; end