# Temperature Advection using finite differences with gridded data

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
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)
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

• can you upload your full code and some sample data somewhere so I can run it locally ? – gansub Mar 22 '19 at 10:14
• @gansub, the full code and sample data can be downloaded from this Dropbox folder: dropbox.com/sh/kcrb08h72jjj3rn/AAAIAUgKVrRrICAQtXhOHS2ta?dl=0 Thank you for your help. – ouranos Mar 22 '19 at 10:43
• Thanks. If you do not mind (now the question is complete). Can you edit that dropbox link into the question ? – gansub Mar 22 '19 at 12:12
• Sure. Done, link is edited into question. – ouranos Mar 22 '19 at 13:03
• The (0,0) position is at 0° lon and 20°N lat (0,20). This corner is not shown in this plot. Actually, all the boudaries are out of the figure domain. – ouranos Mar 25 '19 at 17:15

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.

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.

• I edited the question to remove the $s$ dimention in $T$. It is only a function of time and space $T$[t,y,x]. – ouranos Mar 20 '19 at 19:32
• Three fields are given. T, u and v: temperature, zonal and meridional components of the horizontal wind. That is where I obtained T. I don't understand your comment "...but this seems a bit like a redundant operation: the end result is 𝑇:=𝑇" Could you elaborate more? I am not calculating T. T is given. I am scaling u and v, multiplying each by T, making uT and vT. T is in module, because it contains negative values from anomalies. That is what I called "scaling wind components by temperature". The magnitude of the scaled wind vector is then given by $\sqrt{uT^{2}+vT^{2}}$. – ouranos Mar 20 '19 at 19:39
• Centered finite differences are very commonly used in this field GrADS has been successfuly using it for decades... cola.gmu.edu/grads/gadoc/gradfunchcurl.html You are probably right about the "upwinding". But that should not be the main source of the striking difference between the middle and rightmost figures. I am asking about my assumption that the classical temperature advection should be similar to the magnitude of the temperature-scaled wind. – ouranos Mar 20 '19 at 19:43
• The integration time should not matter, since this is not a prognostic quantity. The question is still valid for one time step. – ouranos Mar 20 '19 at 19:44
• @ouranos: Other than that, if you're still stuck with that problem, I suggest to reduce the size of the problem: See what happens if you would advect only one temperature bump with your scheme along one direction. This would give you a cleaner setup and an easier way to identify what's wrong. – AtmosphericPrisonEscape Mar 22 '19 at 16:05

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.