horizontal_diffusion factor of 2 error

[graig-sutherland]

Hi OpenDrift,

Noticed in OpenDrift that you define random walk velocity to have magnitude np.sqrt(2D/dt), which isn't correct if you are selecting the stochastic component from a normal distribution with standard deviation = 1. The factor of 2 comes from the stochastic differential equation where the variance of stochastic differential is equal to 2dt, but this is not done numerically. I've attached a jupyter notebook which shows the issue. The diffusivity, as calculated from the spread of the particles, will be a factor of 2 greater than the defined value input into OpenDrift (i.e. if 15 m^2/s is input into OpenDrift the true dispersion will be 30 m^2/s). This could be important when trying to input diffusion coefficients from observations into OpenDrift. random_walk_test.pdf

[knutfrode]

Hi Graig, long time - hope you are doing fine!

I will look into the theory again, but do trust that you are right, so that this should be fixed immediately. I guess then the fix should be to simple remove the factor of 2, so that we have:

x_vel = self.elements.moving*np.sqrt(D/dt) * np.random.normal(scale=1, size=self.num_elements_active())
y_vel = self.elements.moving*np.sqrt(D/dt) * np.random.normal(scale=1, size=self.num_elements_active())

https://opendrift.github.io/_modules/opendrift/models/basemodel.html#OpenDriftSimulation.horizontal_diffusion

I believe also @johannesro and @gauteh are interested in this topic.

[graig-sutherland]

Hi Knut-Frode,

It has been a while. All fine here in Canada. Still working from home. I hope this finds you well.

I wanted to flag this as we discovered the same issue in our random walk scheme and a colleague mentioned that you have it coded the same way in OpenDrift. It's due to a disconnect between how the theory is presented and the numerical implementation and as you mention is a trivial thing to fix. This also makes me wonder if this is widespread in various particle tracking codes?

[graig-sutherland]

Hello again,

I found the issue (and it's not OpenDrift). I was adding the diffusion in the two directions together. So the total diffusion (which is what I was calculating) is Kx+Ky, which is where the factor of 2 comes from as Kx=Ky in this case. I think this issue can be closed. Sorry for any confusion.

[knutfrode]

Ok, then all is good. Thank you anyway for warning about the potential issue!

[graig-sutherland]

No problem. If anything this is important for me to remember when comparing diffusion from observations and that from particle tracking.