Temperature Advection using finite differences with gridded data

Question

Advection of a scalar quantity, such as temperature (T), by the horizontal wind, is defined as follows:

$-\textbf{U}\cdot\nabla T$

where $\textbf{U}$ is the horizontal wind vector, $\textbf{U}=(u,v)$, being $u$ and $v$ its zonal and meridional components, respectively (source).

To compute it with gridded data, it is convenient to use finite differences (or similar numerical approaches), as pointed here.

Using this method, the following lines written in Python calculate the advection of temperature (TAdv):

dy=111000 # [m]
lonres=lon[1]-lon[0] # constant
tadv=np.zeros(np.shape(T)); tadv.fill(np.nan)
for t in range(np.shape(T)[0]):
    for x in np.arange(1,len(lon)-1):
        for y in np.arange(1,len(lat)-1):
            dx = abs(111000*np.cos(lat[y]*(2*np.pi/360))*lonres)
            tadv[t,y,x] = -(u[t,y,x]*(T[t,y,x+1]-T[t,y,x-1])/(2*dx) +\
                              v[t,y,x]*(T[t,y+1,x]-T[t,y-1,x])/(2*dy))

I would expect a similar result from scaling the wind components by the temperature, at each grid cell, which is done as follows:

ut  = np.zeros(np.shape(T))
vt  = np.zeros(np.shape(T))
UT  = np.zeros(np.shape(T))
for t in range(np.shape(T)[0]):
    for x in range(len(lon)):
        for y in range(len(lat)):
            ut[t,y,x] = u[t,y,x]*abs(T[t,y,x])
            vt[t,y,x] = v[t,y,x]*abs(T[t,y,x])
            UT[t,y,x] = ((ut[t,y,x]**2 + vt[t,y,x]**2)**0.5)*np.sign(T[t,y,x])

The magnitude of the scaled wind vectors (UT) would be then similar to this "amount of transport of T". It is multiplied by the sign of T (-1 or 1), because T are anomalies (there is T<0), in order to keep its the sign.

However, the results are totally different. What I am missing here? Is it a mistake in this assumption (that both results should look fairly similar)? Or is it a coding problem?

In the figure below, I show the original T, u and v fields (left), and the results of the two above calculations: Tadv with u and v (middle), and the magnitude of the scaled vectors, UT, with the scaled vectors themselves uT and vT (right). TAdv is multiplied by 24*3600, so that the units are K/day.

Here is a link to a Dropbox folder, where the full script and sample data are available for download: https://www.dropbox.com/sh/kcrb08h72jjj3rn/AAAIAUgKVrRrICAQtXhOHS2ta?dl=0

Answer 1

When I look at problems like these I first check to see if there is a well tested and well documented implementation already rather than reinventing the wheel. In this case MetPy temperature advection is a well tested software that does many of the things meteorologists want including calculating finite differences with the right map scale factors. Since you are using Python I would at least look at their implementation first and correct your code even if you are not going to use their code.

Please do this first and validate your data before you attempt any downscaling operation.

And finally AtmosphericPrisonEscape is correct. If you are using lat,lon grid and calculating finite differences you need to include map scale factors. And if you are using global grids the finite differences will not work at the poles. If you are using regional grids that is not a problem. Please include as much detail as possible when you ask a computational question.

Answer 2

So i don't know for how many timesteps you've integrated. And I don't really understand what you have done to obtain the temperature $\rm T[t,s,y,x]$ in your second box, but this seems a bit like a redundant operation: the end result is $T:=T$, so of course your right side looks like the left side, in terms of the temperature colourmap. I don't understand where the different scaling comes from though.

Your first numerical recipe is a whole different beast: Approximating a differential in a first order symmetric fashion as $dT/dx \approx (T[t,y,x+1]-T[t,y,x-1])/(2*dx) $ is known to be numerically unstable. This is why when you want to stick to finite differences, one uses upwinding.

Also your finite differences are cartesian ones, while you project them on a sphere, so to get physically correct results, you need to take the correct differentials in spherical coordinates into account.

Answer 3

I've come to a solution. I have followed @gansub's suggestion to check what there is already there, before reinventing the wheel. I found MetPy's code very (too) well structured, so the code for advection was sparse among several different functions; not so immediate to trace back and put a function together. However, GrADS brings a very nice example on how to easily build such a function in their documentation page for the cdiff function here. So I rewrote their code in Python, and it now works well. I wrote a small package available here, if anyone is interested in the working solution.