Model velocity output validation

[jporobicg]

Hi Everyone, I've not been successfully validating the output from OceanParcel with the hydrodynamic. When I'm comparing the velocity field from the hydrodynamic and the velocity coming from OceanParcel, they don't match. In fact, the particle movement, the velocity (direction) from OP and the velocity from the hydro model are all different (see attached figure). My hydrodynamic model uses a curvilinear c-grid, so I'm using FieldSet.from_nemo. I guess I'm doing something wrong, I got help from the Gitter OP community, but we could not find the problem. I was expecting to have differences but not so big. I'm using the following code to save the U and V that particle 'feels' at each time step.

def AdvectionRK4(particle, fieldset, time):
    """Advection of particles using fourth-order Runge-Kutta integration.
    Function needs to be converted to Kernel object before execution"""
    (u1, v1) = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
    lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt)
    (u2, v2) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat1, lon1]
    lon2, lat2 = (particle.lon + u2*.5*particle.dt, particle.lat + v2*.5*particle.dt)
    (u3, v3) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat2, lon2]
    lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.dt)
    (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3]
    particle.lon += (u1 + 2*u2 + 2*u3 + u4) / 6. * particle.dt
    particle.lat += (v1 + 2*v2 + 2*v3 + v4) / 6. * particle.dt
    particle.U_effective = (u1 + 2*u2 + 2*u3 + u4) / 6.
    particle.V_effective = (v1 + 2*v2 + 2*v3 + v4) / 6.``

Fig1. Comparison of the direction of the velocities of OceanParcels and the hydrodynamic model. a) only de direction; b) direction at particle location Any help will be welcome. Thanks in advance.

[erikvansebille]

Thanks for this question, @jporobicg. We stopped supporting the gitter channel last year; didn't realize there is still activity on it...

Anyways, one reason could be that the velocities that Parcels calculates are 'true' zonal and meridional, while they are in the i- and j-direction in netcdf files. If a grid is rotated, these are not the same. See also the discussion in our parcels v2.0 paper and the Nemo Curvilinear tutorial.

[jporobicg]

Hi @erikvansebille Thanks a lot for the answer and the information; due to the grids' transformations, I will not compare the velocities from OceanParcels and our hydrodynamic model. However, you were right; although the magnitudes of the velocities were similar, the values of U and V between the two models are different, which makes me think that the rotation of the grids caused it (fig 1).

I did compare the particle position and the velocities coming from the hydrodynamic model, and I expected that the direction in which the particles move and the direction of the currents coincide. I understand they won't be the same due to the Runge-Kutta integration, but I didn't expect them to be so different. Am I doing something wrong, or am I not understanding something? (fig2)

[erikvansebille]

Good to see that magnitude in your top plot does agree, that supports my hypothesis that the mismatch between 'U' and 'V' in the GRB1 input and the parcels-output is because of the rotated grid

The same hypothesis could explain the mismatch in the lower figure. You have plotted the quivers on a lat/lon coordinate system, but they should be plotted on the (rotated) grid. See e.g. Figure 6 of the tutorial here for an example how such a grid can be rotated. From looking at the bottom figure, your GRB1 grid is probably oriented approx 45 degrees clockwise?