Identify undulating free surface in vector field

I have the following data set (attached) that cosists of snapshots of a spatial horizontal (data.vx) and vertical (data.vy) components of velocity, along with the respective x and y coordinates (data.x and data.y, respectively) of a wave moving in water. I don't have access to the actual images so this data is all i can use. if you do a quiver plot of the velocity field you can notice that there is a free surface that ondulates and moves from left to right. I'm trying to mark the free surface by fitting a sinusoid to the average of the velocity over a certain depth using mdl = fittype('a*sin(b*x+c)+d','indep','x'); as seen in the following code

close all
[X, Y] = meshgrid(data(1).x,data(1).y);
X=X'; Y=Y'; 
    x = X(:,1);
yesplot = 0;
for i=1:length(data)

    %%
    U =(data(i).vx);
    meanU = nanmean(U(:,25:50)')';

    if yesplot
    plot(x,meanU)
    hold on
    pause(.1)
    end


    mdl = fittype('a*sin(b*x+c)+d','indep','x');

    ind = ~isnan(meanU);
    ind = find(ind);
    ind = ind(5:end-5);
    fst = [.05,2*pi/.1,0,0];
    flo = [0, 2*pi/.5,-2, -.01];
    fup = [.07, 2*pi/.2, 2, .05];
    %options = fitoptions('poly3', 'Normalize', 'on', 'Robust', 'Bisquare')
    fittedmdl = fit(x(ind),meanU(ind),mdl,'start',fst,'lower', flo,'upper',fup);


    a(i) = fittedmdl.a;
    b(i) = fittedmdl.b;
    c(i) = fittedmdl.c;
    d(i) = fittedmdl.d; 
    if yesplot
    plot(fittedmdl)
    hold off
    drawnow
    end
    
    eta = 0.73*a(i)*sin(b(i)*X(:,1)+c(i)) + d(i); % the 0.73 is a guess to scale from velocity (m/s) used for fitting to the surf elevation (m)
end

but it doesn't work most of the time. I also tried binarizing the vector filed and detecting blobs like this

close all
x = X(1,:);
yesplot = 1;
filtwidth = 13;

for i=1:100%length(data)

    vx =(data(i).vx');
    vy =(data(i).vy');

    [vx, vy] = interpolateAndFilter(vx, vy, X', Y', filtwidth);

    idx = abs(vy) < 0.0001 | abs(vx) < 0.0001;
    idx = imfill(~idx, 'holes');
    idx = bwareaopen(idx, 10000);
    idx = ~idx;
    
    yy = flip(data(1).y);
    
    % couldn't get this to work

    [~, topBoundary] = max(flipud(~idx), [], 1);
    topBoundary = yy(topBoundary);
    
    mdl = fittype('a*sin(b*x+c)+d','indep','x');
    fittedmdl = fit(x(5:end-5)',topBoundary(5:end-5)',mdl,'start',[.05,2*pi/.2,0,0],'lower', [0 2*pi/.5 -2 -.01],'upper',[.02 2*pi/.2 2 .05]);

    a(i) = fittedmdl.a;
    b(i) = fittedmdl.b;
    c(i) = fittedmdl.c;
    d(i) = fittedmdl.d;

    if yesplot
    plot(x,topBoundary)
    hold on
    pause(.05)

    end 

    if yesplot
    plot(fittedmdl)
    hold off
    drawnow
    end
    
    eta = a(i)*sin(b(i)*X(:,1)+c(i)) + d(i);
end

but it's even worse. The fitting doesn't need to be perfect but it needs to make sense, that is, the wave phase (parameter c) needs to travel and not jump all over the place, the amplitude needs to grow/shrink slowly, etc.

I'd appreciate any help.