Calculating 3D Frontogenesis from WRF output in Python

Question

I'm currently trying to run a python script to calculate three-dimensional frontogenesis on my WRF output data to perform cross section analysis with. I already have a working two dimensional version of Petterson's Formula in place to compare against, however my three dimensional version is only capturing very minor aspects of what the two dimensional version is.

The formula I'm using to calculate 3D frontogenesis is:

Where the gradient term is:

Here's an example image of what my two dimensional code (925 hPa) produces:

And here's the three dimensional code which is then interpolated down to the same two dimensional (925 hPa) surface:

The fact that some of the readings do at least appear in the same geographical regions on the images show to me that I'm at least on the right track and that I'm probably making a minor error somewhere. From what I've seen, my python code appears correct, at least by my understanding of how the np.gradient function works, here is the code that calculates the frontogenesis:

# Fetch the fields we need
p = getvar(ncFile, "pressure") * 100 # Convert to Pa
z = getvar(ncFile, "z")
ua = getvar(ncFile, "ua")
va = getvar(ncFile, "va")           
theta = getvar(ncFile, "theta")
omega = getvar(ncFile, "omega")

dz = np.gradient(z, axis=0)
dp = np.gradient(p, axis=0)

theta_gradient = np.sqrt((np.gradient(theta, dx, axis=2))**2 + (np.gradient(theta, dy, axis=1))**2 + (np.gradient(theta, axis=0) / dz)**2)

zonal_gradient = (-1 * np.gradient(theta, dx, axis=2)) * ((np.gradient(ua, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dx, axis=2) * np.gradient(theta, dy, axis=1)))
meridional_gradient = (-1 * np.gradient(theta, dy, axis=1)) * ((np.gradient(ua, dy, axis=1) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dy, axis=1) * np.gradient(theta, dy, axis=1)))
vertical_gradient = (-1 * (np.gradient(theta, axis=0) / dp)) * ((np.gradient(omega, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(omega, dy, axis=1) * np.gradient(theta, dy, axis=1)))

F3D = 1.08e9 * (1 / theta_gradient) * (zonal_gradient + meridional_gradient + vertical_gradient)
return F3D

As a reference, the dx and dy terms are obtained directly from the NetCDF file itself (It has attributes DX and DY, both defined as 4000m)

I'm using the wrf-python library to import my data by means of getvar, which is importing the netCDF files. As a reference, the netCDF files use an array structure similar to standard numpy arrays: [bottom-top, south-north, east-west]

So ordering of the axis argument should be correct (z = 0, y = 1, x = 2).

One of my faculty advisers believes the edges of the gradient may be causing some issues with the inside of the calculation, however my understanding is that each point is calculated independently of the edges, and therefore there should be no difference, but this is something I'm not 100% sure about.

Does anyone know why the calculation would produce an incorrect result as shown in the images above?

Answer 1

After some more head bashing, I was able to solve the issue. Goes to show to make sure to consider your full unit analysis!

The term in the gradient of potential temperature was considering the vertical axis as z, while the term in the vertical gradient with respect to omega was considering p. As omega is a vertical unit of pressure, my term selection for the potential temperature gradient was incorrect. Swapping out the derivative from z to p fixed the first half of the problem.

Secondly, when it comes to calculating the derivatives in the vertical direction, you need to consider that numpy.gradient is assuming you have traveled over two data points, and thus it divides the result by two incorrectly, so I created a function wrapper to handle the inner and outer points of the analysis:

def calc_center_difference(A, ax):
    gradient = np.gradient(A, axis=ax) 
    gradient *= 2.0
    if ax == 0:
        gradient[0,:,:] /= 2.0
        gradient[-1,:,:] /= 2.0   
    elif ax == 1:
        gradient[:,0,:] /= 2.0
        gradient[:,-1,:] /= 2.0
    elif ax == 2:
        gradient[:,:,0] /= 2.0
        gradient[:,:,-1] /= 2.0
    else:
        return ValueError("Invalid axis passed to calc_center_difference")
    return gradient

Afterwards, I swapped the instances of derivative terms that considered either dz or dp with the wrapper function instance to make sure the terms were being correctly calculated, and voila! Everything works!

Here's the corrected form of my earlier function that can calculate the full 3D frontogenesis in python:

# Input netcdf
# - [bottom_top, north_south, west_east]
def three_d_fronto(ncFile):
    # Fetch the fields we need
    p = to_np(getvar(ncFile, "pressure") * 100)
    z = to_np(getvar(ncFile, "z"))
    ua = to_np(getvar(ncFile, "ua"))
    va = to_np(getvar(ncFile, "va"))    
    theta = to_np(getvar(ncFile, "theta"))
    omega = to_np(getvar(ncFile, "omega"))

    dz = calc_center_difference(z, 0)
    dp = calc_center_difference(p, 0)

    theta_gradient = np.sqrt((np.gradient(theta, dx, axis=2))**2 + (np.gradient(theta, dy, axis=1))**2 + (calc_center_difference(theta, 0) / dp)**2)
    zonal_gradient = (-1 * np.gradient(theta, dx, axis=2)) * ((np.gradient(ua, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dx, axis=2) * np.gradient(theta, dy, axis=1)))
    meridional_gradient = (-1 * np.gradient(theta, dy, axis=1)) * ((np.gradient(ua, dy, axis=1) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dy, axis=1) * np.gradient(theta, dy, axis=1)))
    vertical_gradient = (-1 * (calc_center_difference(theta, 0) / dp)) * ((np.gradient(omega, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(omega, dy, axis=1) * np.gradient(theta, dy, axis=1)))

    F3D = 1.08e9 * (1 / theta_gradient) * (zonal_gradient + meridional_gradient + vertical_gradient)
    return F3D

Just as a reminder, you need the wrf-python library installed to use this function.