Passive advection and comparison with parcels

[margauxf]

1. In an effort to model the trajectories of drogues (SVP-type), I am writing my own model for 2-D advection of 3-D objects. I have followed openberg, opendrift and shipdrift, if you have any other suggestion, I'd love to hear it.

2. I have been comparing trajectories with the output from oceanparcels but keep running into discrepancies in the results, even with same initial conditions (lon, lat, time). Because the parcels package only advects passive trajectories, I tried to advect OpenDrift with a wind drift factor of 0:

   # Second run, for comparison
   o2 = OceanDrift(loglevel=20)  # Set loglevel to 0 for debug information
   o2.add_reader([reader_current])
   o2.seed_elements(lon = -65, lat = 40., time = t_i, wind_drift_factor=0.0)
   o2.set_config('drift:scheme', 'runge-kutta4')

And from the model template, I created a new model that has nothing but the following:

self.update_positions(self.environment.x_sea_water_velocity, self.environment.y_sea_water_velocity)
self.deactivate_elements(self.environment.land_binary_mask==1)

yet I still see discrepancies (see figure).

How can I ensure I am modeling purely passive particles? Do you know the reasons for these discrepancies, and the different results from parcels? I ensured I used RK4 for all simulations.

Thanks for a great toolbox, the embedded HYCOM readers are such time savers.

[knutfrode]

1. It should not be necessary to write a new module for this purpose, as this can be achieved with OceanDrift. As you might have seen, OceanDrift3D was merged into OceanDrift yesterday. So OceanDrift is now a 3D model, but any vertical movement may be switched off with the configuration mechanism. I now also made a dedicated method set_2d() as a shortcut convenience for this. A new example of this method should appear here within an hour or two. You should also remember to seed your SVP-particles at correct depth, z=-15.

2. It is very interesting that you have compared OpenDrift to OceanParcels. In fact, I would say that the comparison you show is rather good. The drift problem is highly unstable, and any tiny error (numerical or methodological) would quickly grow. Although conditions seem to be identical, there are several choices inside each model which have some numerical impact: horizontal interpolation, vertical interpolation, temporal interpolation, time step, great circle calculation, reprojection, vector rotation.... Also note the interesting paper of Nordam et al. which says that the integrator (Euler, RK2, RK4) makes little difference if the interpolator is linear. You may also have seen this example script illustrating error dependency on integrator and calculation time step. It is of course not an excuse to be sloppy with accuracy, but as found in my own study, the difference between a real and simulated drifter is much larger (500m offset per hour) than the difference you show. Thus most of the error is in the ocean model (or probably its (limited) data assimilation), and not in the drift model.

[margauxf]

Thank you for getting back to me promptly!

1. I will have a look at the new release of OceanDrift3D. That looks very promising. One of the reasons why I was looking into creating a new module was to input specific values of drag/area for certain objects, especially at depth. This is why I was looking at openberg (the area coefficient varies with depth) and shipdrift (the form drag is incorporated here). I am uncertain that wind drift values are sufficient proxies to model 2-D advection of objects with certain form drag. Perhaps you have some insight into this.

2. Thank you! I am sorry for not being clearer, the picture was NOT comparing OpenDrift to Ocean Parcels, but it was comparing a totally passive advection (with a new module) to setting wind drift = 0. This is the passive particle module I called:

import numpy as np
from opendrift.elements import LagrangianArray
from opendrift.models.basemodel import OpenDriftSimulation

class TemplateElementType(LagrangianArray):
    variables = LagrangianArray.add_variables([
        ('wind_drift_factor', {'dtype': np.float32,
                               'unit': '%',
                               'default': 0.00}),
        ('terminal_velocity', {'dtype': np.float32,
                               'units': 'm/s',
                               'default': 0.})])

class PassiveTest(OpenDriftSimulation):
    ElementType = TemplateElementType

    required_variables = ['x_sea_water_velocity', 'y_sea_water_velocity',
                          'x_wind', 'y_wind', 'land_binary_mask']

    fallback_values = {'x_sea_water_velocity': 0,
                       'y_sea_water_velocity': 0,
                       'x_wind': 0, 'y_wind': 0}

    def __init__(self, *args, **kwargs):
        super(PassiveTest, self).__init__(*args, **kwargs)

    def update(self):
        self.update_positions(self.environment.x_sea_water_velocity,
            self.environment.y_sea_water_velocity)
        self.deactivate_elements(self.environment.land_binary_mask==1)

and I compared the two:

# Simulation 1
o = PassiveTest(loglevel=20) 
o.add_reader([reader_current])
o.seed_elements(lon = -65, lat = 40., time=t_i)
o.set_config('drift:scheme', 'runge-kutta4')
o.run(steps=31, time_step=3600, outfile = 'opendrift_output_drifters_passivetest.nc')

# Simulation 2
o2 = OceanDrift(loglevel=20)  # Set loglevel to 0 for debug information
o2.add_reader([reader_current])
o2.seed_elements(lon = -65, lat = 40., time = t_i, wind_drift_factor=0.0)
o2.set_config('drift:scheme', 'runge-kutta4')
o2.run(steps=31, time_step=3600, outfile = 'opendrift_output_drifters_0.0wind_RK4.nc')

While they are similar, I was wondering about the discrepancy. Is it because OpenDrift incorporates stokes drift? With -what I thought - were the exact same parameters, I was expecting the two trajectories to be identical.

I agree with you that if the difference between Ocean Parcels and OpenDrift was this small, it would be very good. Please see the newly uploaded picture for comparison. I verified the timestamps. It looks as though there is some drag to the OpenDrift trajectories in comparison to the Parcels trajectory.

as found in my own study, the difference between a real and simulated drifter is much larger (500m offset per hour) than the difference you show. Thus most of the error is in the ocean model (or probably its (limited) data assimilation), and not in the drift model.

I absolutely agree with you and I have read your paper. It is my view, in my work, that perfecting drag/added mass coefficients is not worth the time for advection studies on the scale of months, as very few models a) have high enough resolution over such large domains and b) have high enough accuracies. My goal was to rather have an understanding of the spread of trajectories (and ultimately output a spectrum of dispersion) to understand the role of drag.

Thank you again for your time.

[knutfrode]

Thank you for getting back to me promptly!

1. I will have a look at the new release of OceanDrift3D. That looks very promising. One of the reasons why I was looking into creating a new module was to input specific values of drag/area for certain objects, especially at depth. This is why I was looking at openberg (the area coefficient varies with depth) and shipdrift (the form drag is incorporated here). I am uncertain that wind drift values are sufficient proxies to model 2-D advection of objects with certain form drag. Perhaps you have some insight into this.

I have not done any development on OpenBerg or ShipDrift models myself, and in all other models it is assumed that all objects follow current 100%, plus Stokes (also 100%) and a percentage of the wind. The wind_drift_factor is anyway empirical, so I guess the form drag is implicit there, as also in the Leeway coefficients. But perhaps using an explicit form drag would give a percentage which depends on the wind speed? I don't have much insight into the impact of a varying form drag, so any contributions to the codebase from you or others on this would be welcome! Eventually it would be good to incorporate this into the generic OceanDrift model, but for the development it might perhaps be easier to start with a minimalistic model as you have already done.

2. Thank you! I am sorry for not being clearer, the picture was NOT comparing OpenDrift to Ocean Parcels, but it was comparing a totally passive advection (with a new module) to setting wind drift = 0. This is the passive particle module I called:

import numpy as np
from opendrift.elements import LagrangianArray
from opendrift.models.basemodel import OpenDriftSimulation

class TemplateElementType(LagrangianArray):
    variables = LagrangianArray.add_variables([
        ('wind_drift_factor', {'dtype': np.float32,
                               'unit': '%',
                               'default': 0.00}),
        ('terminal_velocity', {'dtype': np.float32,
                               'units': 'm/s',
                               'default': 0.})])

class PassiveTest(OpenDriftSimulation):
    ElementType = TemplateElementType

    required_variables = ['x_sea_water_velocity', 'y_sea_water_velocity',
                          'x_wind', 'y_wind', 'land_binary_mask']

    fallback_values = {'x_sea_water_velocity': 0,
                       'y_sea_water_velocity': 0,
                       'x_wind': 0, 'y_wind': 0}

    def __init__(self, *args, **kwargs):
        super(PassiveTest, self).__init__(*args, **kwargs)

    def update(self):
        self.update_positions(self.environment.x_sea_water_velocity,
            self.environment.y_sea_water_velocity)
        self.deactivate_elements(self.environment.land_binary_mask==1)

and I compared the two:

# Simulation 1
o = PassiveTest(loglevel=20) 
o.add_reader([reader_current])
o.seed_elements(lon = -65, lat = 40., time=t_i)
o.set_config('drift:scheme', 'runge-kutta4')
o.run(steps=31, time_step=3600, outfile = 'opendrift_output_drifters_passivetest.nc')

# Simulation 2
o2 = OceanDrift(loglevel=20)  # Set loglevel to 0 for debug information
o2.add_reader([reader_current])
o2.seed_elements(lon = -65, lat = 40., time = t_i, wind_drift_factor=0.0)
o2.set_config('drift:scheme', 'runge-kutta4')
o2.run(steps=31, time_step=3600, outfile = 'opendrift_output_drifters_0.0wind_RK4.nc')

While they are similar, I was wondering about the discrepancy. Is it because OpenDrift incorporates stokes drift? With -what I thought - were the exact same parameters, I was expecting the two trajectories to be identical.

Ok, it was too good to be true with such small difference among two different model frameworks (OceanParcels and OpenDrift) :-) I think I realise now what makes the difference between your minimalistic model and OceanDrift:

• one thing might be that OceanDrift might still add Stokes drift, if there is wind in your forcing input. You could make sure to switch off Stokes drift by o.set_config('drift:stokes_drift', False) You could also check the log from your run (loglevel=0) to see everything that goes on.
• but most importantly is that your module has used the low level method update_positions for the advection. This method only gives a pure Euler propagation. I guess some parts of the documentation is slightly outdated, but instead you should use the higher level advection methods found in physics_methodsmodule, e.g. for current: https://github.com/OpenDrift/opendrift/blob/master/opendrift/models/physics_methods.py#L277 Note also that this, and several other of the advection methods, have an input factor (default 1) to scale the movement. So you may e.g. use a factor of 0.8 to move with 80% of the current. Presently this is only used for oil in ice, which oil moves with a fraction of current and ice velocities, etc.

I agree with you that if the difference between Ocean Parcels and OpenDrift was this small, it would be very good. Please see the newly uploaded picture for comparison. I verified the timestamps. It looks as though there is some drag to the OpenDrift trajectories in comparison to the Parcels trajectory.

Ok, this difference with Parcels is more substantial, and should be understood. I agree it looks like there is an additional drag to the OceanDrift results, which could be windage or Stokesdrift, but this is anyway not included in your PassiveTracer model. From your code samples, the OpenDrift simulations are truly on the very surface (z=0), and this is presumably also the case for OceanParcels? (typically the internal loops are more prominent at depth, but smeared out by windage closer to the surface). Could you eventually repeat the simulation with z=-15 for the run with OceanDrift?

[margauxf]

Thank you very much! This is very helpful. I will run the tests you suggested and get back to you. I think you are right about the Stokes drift.

I don't have much insight into the impact of a varying form drag, so any contributions to the codebase from you or others on this would be welcome!

Fun! I will share my module once I'm done with it then :-)

Eventually it would be good to incorporate this into the generic OceanDrift model, but for the development it might perhaps be easier to start with a minimalistic model as you have already done.

Yes, I wanted to make sure I understood the advection before having:

• Stokes drift, windage impacting the surface
• subsurface currents impacting the water column ~[-2, -15]m

[knutfrode]

I managed to reproduce your OpenDrift simulations using HYCOM data from http://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0/uv3z And in searching for the difference with your results using Parcels, I found and fixed a severe bug in the netCDF-reader, which gave wrong results west of the 0-meridian for the HYCOM dataset (longitudes 0-360).

After this fix, the OpenDrift simulation now coincides very well with what you got with Parcels. So thank you very much for your help in spotting this one.

Btw, in the meantime I also updated the documentation in the model_template, which was rather outdated. It should still be improved, but at least it is better than before.

[margauxf]

Fantastic news! Indeed, the agreement is remarkable. I can now keep building a new model template. Thank you so much for your help!

[margauxf]

Hi again,

Based on the drifter example I tried reiterating and comparing the results by reading HYCOM data:

# Simulation 1
o = OceanDrift(loglevel=20) 
o.add_readers_from_list(['https://thredds.met.no/thredds/dodsC/sea/norkyst800m/1h/aggregate_be'])

# Simulation 2
o2 = OceanDrift(loglevel=20)
data_url = "http://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0/uv3z"
reader_current = reader_netCDF_CF_generic.Reader(data_url)
o2.add_reader([reader_current])

When seeding trajectories, however, it seems as though the different wind drift factors simply translate the trajectories as opposed to impacting their dynamics.

This is not observed with the NorKyst data:

I tried downloading the .nc files from the HYCOM server to check it wasn't the reader but the y spacing isn't uniform.

I don't think it has to do with the resolution of the HYCOM data (0.083 degrees as opposed to NorKyst's 800m) because I am not seeing this issue with OceanParcels. Granted, I have not incorporated drag nor wind drift in the Parcels trajectories.

Thanks in advance!

[knutfrode]

I guess the difference here is that the NorKyst dataset includes both current and wind, whereas the HYCOM dataset you are using includes only currents. Thus a fallback value of 0 is then used for wind. If you change loglevel to 0, you will see messages about this in the log.

You could try to add also this NCEP dataset along with HYCOM to get winds as well: https://pae-paha.pacioos.hawaii.edu/thredds/dodsC/ncep_global/NCEP_Global_Atmospheric_Model_best.ncd

The small translation you see is probably due to the initial random perturbations from a seeding radius of 3000m in that example.