Comparison of divergence algorithms

TL;DR: Be careful when calculating spatial derivatives. Remember to include cosine weighting, convert the unit into radians, and use the full wind instead of the ageostrophic wind only.

I have seen two common algorithms for calculating divergence from wind.
One is from Windspharm:

1
2
3
4

##Suppose u, v are variables for u wind and v wind and have only two dimensions (latitude and longitude)
w = VectorWind(u, v)
div3 = w.divergence()
div3 #This is the variable for divergence

another is from Metpy.

1
2
3

## Again, Suppose u, v are variables for u wind and v wind and have only two dimensions (latitude and longitude)
div4 = mpcalc.divergence(u, v)
div4 #This is the variable for divergence

I also tried to calculate it myself.

1
2
3
4
5
6
7
8
9
10

a = 6378000
cphi = np.cos(data.lat*np.pi/180)
dx = np.gradient(data.lon, axis=0, edge_order=2)/180*np.pi
dx = np.outer(cphi.values, dx)
dx = xr.DataArray(dx*a, dims = ['lat','lon'],
                    coords=dict(lon = data.lon, lat = data.lat))
dy = np.pi/180*a
div_u = ((u*cphi).differentiate('lon', edge_order=2)/(dx*cphi*4)) #/4 is important
div_v = ((v*cphi).differentiate('lat', edge_order=2)/(dy*cphi))
div2 = div_u + div_v

(It worth to mention why /4 is important there. We do that since u.differentiate will have the results with the unit in 1/deg, but out dx is calculated to have the unit in meters/0.25 degrees, so /4 is required to complete the conversion)

However, the mothod above is quite stupid and redundant. The better way is shown below:

1
2
3
4
5

a = 6378000
cphi = np.cos(sample.lat*np.pi/180)
div_u = (test_u).differentiate('lon', edge_order=2)
div_v = (test_v*cphi).differentiate('lat', edge_order=2)
div2 = (div_u + div_v)/cphi/(a*units('meter'))*(180/np.pi)

Here is the reason:
A cautious way to calculate spatial derivatives would be from the original definition from the Cartesian coordinate:

$$Div(V)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}$$

Where:

δx=a∗cosϕ∗δλ,
δy=a∗δϕ

a: Radius of Earth,
ϕ: Latitude,
λ: Longitude.

However, it is also required to correct the result by multiplying them with the cosines, and the equation becomes.

$$Div(V)=\frac{1}{cos(\phi)}\left(\frac{\partial (cos(\phi)u)}{acos(\phi) \partial \lambda}+\frac{\partial (cos(\phi)v)}{a\partial \phi}\right)$$

And you can derivative the equation into:

$$Div(V)=\frac{1}{a*cos(\phi)}\left( cos(\phi)\frac{\partial u}{cos(\phi) \partial \lambda} + u\frac{\partial cos(\phi)}{cos(\phi) \partial \lambda} +\frac{\partial (cos(\phi)v)}{\partial \phi}\right)$$

After simple simplifications, it becomes the classic equation for computing divergence on the spherical coordinate.

$$Div(V)=\frac{1}{a*cos(\phi)}\left(\frac{\partial u}{\partial \lambda}+\frac{\partial (cos(\phi)v)}{\partial \phi}\right)$$

I also tried other ways, one way is to calculate without adding the cosines, which means instead of calculating by:

1

div4 = mpcalc.divergence(u, v)

The calculation is rather done as:

1

div5 = mpcalc.first_derivative(u, axis='lon') + mpcalc.first_derivative(v, axis='lat')

The other is to compute the divergence by ageostrophic wind. However, ERA5 does not have native geostrophic wind variables, so I have to calculate from the geopotential height.

1
2
3
4

ug, vg = mpcalc.geostrophic_wind(z)
div_u = u - ug
div_v = v - vg
div6 = mpcalc.divergence(div_u, div_v)

For the sake of making the results comparable to each other and have a good referee, I, again, use the ERA5 data from the data store, and I am treating the divergence data from the dataset as the referee. I guess I have faith on the ERA5 maintainers… By the way, the divergence variable from the dataset should be used as dataset.d (instead of dataset.divergence, and I wonder if this is as a result of a laziness)…

Here are some results:

It appears like they are all similar in the midlatitudes, but the ageostrophic version seems to have extreme values in the tropics. Windspharm and the method without cosine weighting, respectively, have issues near the poles.
To make the difference more obvious, I made the plots below to compare the results from different methods. Please note that the upper-left one is the reference plot, so it is not a result of subtraction.

The issue in the ageostrophic version should be from calculation of ageostrophic wind in the tropics. This is not hard to understand but unintuitive to me. In short, the low Coriolis wind in the tropics causes wind flow not under geostrophic balance. The low geopotential gradient also leads to poor evaluation accuracy of geostrophic wind. These two issues combined causes the extreme values in the tropics.

Also, the computation time for each method is displayed at the upper right of each figure. Windspharm has not only the worst performance, among all the regularly working methods, but also the slowest computation speed. The native method is the quickest, and the metpy function spends 10x time than the native method. Since both methods are similarly accurate, the suggestion is to calculate natively.

In terms of conclusion? I cannot think about any, but we probably should recommend people to use the ERA5 data directly if possible. I probably should write more about the algorithm behind each method, but I don’t have that much energy to do that at this moment. Let here be the end then.