Check factor of 2 in Waves/SWIFT_Stokes.m

[jacobrdavis]

Code was revised in Mar 2021 to include factor of 2 at lines 176-183, however comparison to other published equations suggest this factor of 2 is not necessary. For example, see Equation 5 in this paper (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020MS002172) or Equation 9.9 in Fabrice's textbook (https://github.com/ardhuin/waves_in_geosciences/blob/main/part1and2.pdf).

Confusion may stem from the "2" in both the numerator and denominator of the integrand:

I.e., to match equations, a factor of "2" needs to be added to the denominator in lines 164-172.

SWIFT-codes/Waves/SWIFT_Stokes.m

Lines 156 to 184 in f07c03d

	% Loop over frequency
	for jj=1:length(frequency)
	om = 2*pi.*frequency(jj);
	% Compute Stokes drift and Wave Bias for this frequency
	if strcmp(depthflag,'intermediate')==1
	k = wavenumber(frequency(jj), h);
	Stokes_jj_0 = om*k*cosh(2*k*h)./(sinh(k*h).^2)*energy(jj);
	Stokes_jj = om*k*cosh(2*k*(z+h))./(sinh(k*h).^2)*energy(jj);
	WaveBias_jj_0 = om*k*cosh(2*k*h)./(sinh(k*h).^2)*energy(jj);
	WaveBias_jj = om*k*cosh(k*z+2*k*h)./(sinh(k*h).^2)*energy(jj);
	else
	Stokes_jj_0 = om*k*energy(jj);
	Stokes_jj = om*k*energy(jj)*exp(2*k*z);
	WaveBias_jj_0 = om*k*energy(jj);
	WaveBias_jj = om*k*energy(jj)*exp(k*z);
	end
	% Sum component-wise with contributions from other frequencies ** added factors of 2 in Mar 2021
	Stokes_spectral_east_0 = Stokes_spectral_east_0+2*Stokes_jj_0.*(-a1(jj))*df;
	Stokes_spectral_north_0 = Stokes_spectral_north_0+2*Stokes_jj_0.*(-b1(jj))*df;
	Stokes_spectral_east = Stokes_spectral_east+2*Stokes_jj.*(-a1(jj))*df;
	Stokes_spectral_north = Stokes_spectral_north+2*Stokes_jj.*(-b1(jj))*df;
	WaveBias_spectral_east = WaveBias_spectral_east+2*WaveBias_jj.*(-a1(jj))*df;
	WaveBias_spectral_north = WaveBias_spectral_north+2*WaveBias_jj.*(-b1(jj))*df;
	WaveBias_spectral_east_0 = WaveBias_spectral_east_0+2*WaveBias_jj_0.*(-a1(jj))*df;
	WaveBias_spectral_north_0 = WaveBias_spectral_north_0+2*WaveBias_jj_0.*(-b1(jj))*df;
	end

Here is a MATLAB test case using a portion of the code above (please double check). Using my own code (not shown), I calculate a value that is half of the value output by this code: a Stokes magnitude of 0.0219 compared to 0.0407 (differences may be due to dispersion calculation/numerical integration).

% test comparison of .m and .py stokes drift codes

energy = [0.00174398, 0.00273885, 0.00324009, 0.00126474, 0.004198  , ...
          0.00838162, 0.00286359, 0.00393084, 0.00944725, 0.01648753, ...
          0.02518561, 0.01375298, 0.01577215, 0.03096342, 0.05275692, ...
          0.09791758, 0.24107521, 0.31352153, 0.31425553, 0.27684291, ...
          0.18990934, 0.14168773, 0.10859684, 0.1049892 , 0.09133319, ...
          0.0688494 , 0.06641062, 0.05996119, 0.04833707, 0.04299529, ...
          0.04031001, 0.0391111 , 0.04241879, 0.03477004, 0.03336525, ...
          0.02650127, 0.02101617, 0.01656881];

frequency = [0.0293 , 0.03906, 0.04883, 0.05859, 0.06836, 0.07813, 0.08789, ...
             0.09766, 0.10742, 0.11719, 0.12695, 0.13672, 0.14648, 0.15625, ...
             0.16602, 0.17578, 0.18555, 0.19531, 0.20508, 0.21484, 0.22461, ...
             0.23438, 0.24414, 0.25391, 0.26367, 0.27344, 0.2832 , 0.29297, ...
             0.30273, 0.3125 , 0.32227, 0.33203, 0.35156, 0.38086, 0.41016, ...
             0.43945, 0.46875, 0.49805];

z = 0;

h = 22; 

a1 = [0.012708,  0.009775, -0.079179, -0.074291,  0.201369,  0.2913  , ...
      0.159335,  0.27175 ,  0.236559,  0.301075,  0.30303 ,  0.243402, ...
      0.215054,  0.217986,  0.231672,  0.226784,  0.230694,  0.233627, ...
      0.222874,  0.166178,  0.1652  ,  0.192571,  0.184751,  0.115347, ...
      0.13783 ,  0.121212,  0.105572,  0.086022,  0.049853,  0.043988, ...
      0.000978, -0.052786, -0.081134, -0.259042, -0.394917, -0.404692, ...
     -0.478006, -0.502444];

b1 = [0.029326, -0.029325, -0.030303,  0.030303,  0.211144,  0.399805, ...
      0.274682,  0.317693,  0.457478,  0.600196,  0.753666,  0.690127, ...
      0.674487,  0.800587,  0.855327,  0.879765,  0.903226,  0.892473, ...
      0.889541,  0.894428,  0.874878,  0.861193,  0.845552,  0.887586, ...
      0.865103,  0.835777,  0.831867,  0.831867,  0.822092,  0.823069, ...
      0.812317,  0.800587,  0.756598,  0.653959,  0.565005,  0.44868 , ...
      0.397849,  0.370479];


[stokes_drift_east, stokes_drift_north] = stokes_drift(energy, frequency, z, h, a1, b1)

function [Stokes_spectral_east, Stokes_spectral_north] = stokes_drift(energy, frequency, z, h, a1, b1)
    % Compute Spectral Stokes and Wave Bias Estimates (using Fourier moments)
    
    depthflag = 'intermediate';

    % Initialize with zero value prior to summing component-wise over freq
    Stokes_spectral_east_0 = 0; % Surface value
    Stokes_spectral_north_0 = 0; % Surface value
    Stokes_spectral_east = zeros(size(z)); % At mean z of Signature bins
    Stokes_spectral_north = zeros(size(z)); % At mean z of Signature bins
    WaveBias_spectral_east_0 = 0; % Surface value
    WaveBias_spectral_north_0 = 0; % Surface value
    WaveBias_spectral_east = zeros(size(z)); % At mean z of Signature bins
    WaveBias_spectral_north = zeros(size(z)); % At mean z of Signature bins

    % Frequency resolution (*assumes constant*)
    df = mean(abs(diff(frequency)));
    
    % Loop over frequency
    for jj=1:length(frequency)
        
        om = 2*pi.*frequency(jj);
        
        % Compute Stokes drift and Wave Bias for this frequency
        if strcmp(depthflag,'intermediate')==1
            k = wavenumber(frequency(jj), h);
            Stokes_jj_0 = om*k*cosh(2*k*h)./(sinh(k*h).^2)*energy(jj);
            Stokes_jj = om*k*cosh(2*k*(z+h))./(sinh(k*h).^2)*energy(jj);
            WaveBias_jj_0 = om*k*cosh(2*k*h)./(sinh(k*h).^2)*energy(jj);        
            WaveBias_jj = om*k*cosh(k*z+2*k*h)./(sinh(k*h).^2)*energy(jj);
        else
            Stokes_jj_0 = om*k*energy(jj);
            Stokes_jj = om*k*energy(jj)*exp(2*k*z);
            WaveBias_jj_0 = om*k*energy(jj);        
            WaveBias_jj = om*k*energy(jj)*exp(k*z);
        end
        
        % Sum component-wise with contributions from other frequencies ** added factors of 2 in Mar 2021
        Stokes_spectral_east_0 = Stokes_spectral_east_0+ 2*Stokes_jj_0.*(-a1(jj))*df; 
        Stokes_spectral_north_0 = Stokes_spectral_north_0 + 2*Stokes_jj_0.*(-b1(jj))*df; 
        Stokes_spectral_east = Stokes_spectral_east + 2*Stokes_jj.*(-a1(jj))*df; 
        Stokes_spectral_north = Stokes_spectral_north + 2*Stokes_jj.*(-b1(jj))*df; 
        WaveBias_spectral_east = WaveBias_spectral_east + 2*WaveBias_jj.*(-a1(jj))*df;
        WaveBias_spectral_north = WaveBias_spectral_north + 2*WaveBias_jj.*(-b1(jj))*df;
        WaveBias_spectral_east_0 = WaveBias_spectral_east_0 + 2*WaveBias_jj_0.*(-a1(jj))*df;
        WaveBias_spectral_north_0 = WaveBias_spectral_north_0 + 2*WaveBias_jj_0.*(-b1(jj))*df;    
    end


end