Q_Vector returns wrong units

[sodoesaburningbus]

What went wrong?

mpcalc.q_vector returns values with units of m2/kg/s. Appropriate units should be kg/m2/s2.

Operating System

Windows

Version

1.6.2

Python Version

3.12

Code to Reproduce

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
import metpy.constants as mpconstants
from metpy.units import units
import numpy as np
import xarray as xr
from datetime import datetime, timedelta, time 

e = datetime(2023, 5, 1, 12)
ds = xr.open_dataset(f'https://www.ncei.noaa.gov/thredds/dodsC/model-gfs-g4-anl-files/'
                     f'{date:%Y%m}/{date:%Y%m%d}/gfs_4_{date:%Y%m%d_%H}00_000.grb2')

# Sunset geographic extent of data to CONUS to limit how much is downloaded
lon_slice = slice(-125, -67)  
lat_slice = slice(24, 50)     

lats = ds['lat'].values  # Extract latitude
lons = ds['lon'].values  # Extract longitude

# Create meshgrid from 1D lat and lon arrays
lons_grid, lats_grid = np.meshgrid(lons, lats)

# Compute grid spacings (dx, dy) using np.gradient
dy, dx = np.gradient(lats_grid, axis=0), np.gradient(lons_grid, axis=1)

# Grab data and smooth using a nine-point filter applied 50 times to grab the synoptic signal
level = 700 * units.hPa
hght_700 = ds['Geopotential_height_isobaric'].metpy.sel(vertical=level)
tmpk_700 = ds['Temperature_isobaric'].metpy.sel(vertical=level)
uwnd_700 = ds['u-component_of_wind_isobaric'].metpy.sel(vertical=level)
vwnd_700 = ds['v-component_of_wind_isobaric'].metpy.sel(vertical=level)
w_700 = ds['Vertical_velocity_geometric_isobaric'].metpy.sel(vertical=level)

# Compute the Q-vector components
uqvect, vqvect = mpcalc.q_vector(uwnd_700, vwnd_700, tmpk_700, level)
print(uqvect)
print(uqvect.to_base_units())
# Compute the divergence of the Q-vectors calculated above (term A of the RHS according to Hoskins et al. (1978)
q_div = -2 * mpcalc.divergence(uqvect, vqvect)

Errors, Traceback, and Logs

<xarray.DataArray (time: 1, lat: 361, lon: 720)> Size: 2MB
<Quantity([[[ 2.11213126e-08 -1.56944484e-10 -1.30758030e-10 ... -9.80448728e-11
   -1.17662229e-10 -1.37315293e-10]
  [ 1.75598535e-10  1.84654699e-10  5.25826480e-09 ...  2.46977001e-10
    2.28100511e-10  2.36409042e-10]
  [ 2.77543024e-12  1.05006759e-12  1.03417809e-12 ...  1.06336328e-12
    1.83656116e-12  1.83360512e-12]
  ...
  [ 5.59318956e-12  8.84133153e-13  2.63383217e-13 ...  5.14350542e-13
    9.62793471e-13  1.80597172e-12]
  [-1.91209715e-09  4.78330514e-10  5.06398520e-10 ... -8.15684653e-10
   -7.88212573e-10 -7.62894144e-10]
  [ 2.62437847e-10 -2.68241094e-10 -2.87836800e-10 ... -2.29051163e-10
    2.82836771e-13  3.12410213e-13]]], 'meter ** 2 / kilogram / second')>

Units at end of array are incorrect.

[kgoebber]

@sodoesaburningbus thanks for submitting this bug report. The error is in the static stability parameter that is the default value doesn't contain any units, even though the documentation states that the parameter should be a value with units, the default value violates that documentation!

[dopplershift]

@kgoebber Good catch and thanks @sodoesaburningbus for the report. I'll note that the reason this isn't currently erroring is because the unit-checking decorator isn't being given a dimensionality for static stability. I'll knock out a quick patch.

[dopplershift]

The current test test_q_vector_without_static_stability is actually specifically looking for these units. I think it was an oversight when we merged it 6 years ago and should have made the default correct on units. But that explains how this has passed as is.

[jthielen]

I needed to check back with #661 and the textbooks to confirm this, but from what I recalled the unit difference was intended behavior, as it was accounting for differing definitions of Q-Vector (between Holton and Bluestein...Holton doesn't include the 1/sigma factor). Of course, over-6-year-ago me didn't clarify this in the docstring, so the mismatch in expectations is certainly understandable!

Though, @sodoesaburningbus, I did want to check: you mentioned expected units as kg/m2/s2. Our tests for the full version (accounting for static stability) have kg/m2/s3. Were you indeed expecting kg/m2/s2? If so, could you point us to a reference with that information, since our prior sources led (at least me) to believe it should be kg/m2/s3

[dopplershift]

I have #3690, but I do think we need to think hard about the default for the static stability parameter to q_vector. I don't love changing units, but I'm also not convinced that passing in a value of 1 with proper units is a good idea either. Unfortunately, changing the value starts to make this a bigger impact.

[kgoebber]

I now also took a look at Holton (4th ed.) and what is done there is keeping the static stability parameter on the operator side of the equation, which is different than what was done in the traditional QG-omega equation. Why he didn't divide that through in the same fashion to define Q is left unknown. In the Holton and Hakim (5th ed.), Q is derived in a substantially different manner (as is done in other parts of the text as well) from what had been the case through the first four editions. The equation is geared more toward a constant height as opposed to a constant pressure derivation.

Lackmann (2011) follows Bluestein, whereas Martin (2005) follows Holton (4th ed). So clear as mud. I think the leaving of the static stability parameter on the operator side is more done of convenience in writing the equation. In the end, static stability is one of the largest factors in modulating rising motion that is not either vorticity advection, temperature advection, or diabatic processes. As a result, this would be why I would advocate for using static stability with units in our function as opposed to not using it. Additionally, since we publish the Bluestein equation, which includes the static stability parameter, it seems logical that we would produce the calculation with the appropriate units for that particular equation.

[sodoesaburningbus]

Thank you for the rapid response. It was actually one of my students who caught the bug.

[sgdecker]

It might also be worth checking to see how GEMPAK handles (or doesn't handle) the static stability parameter, if we are looking for feature parity. Or perhaps an optional argument to choose which flavor of the Q-vector to compute.

[sgdecker]

I've looked into this in more detail. First, GEMPAK has two Q-vector functions, QVEC, which doesn't include the static stability parameter, and QVCL, which does. Second, in the original Hoskins et al. (1978) paper that introduced Q-vectors (which uses a similar vertical coordinate as Holton/Hakim 5th ed), the static stability parameter ($N^2$ in this coordinate) is kept on the left-hand side (their Eq. 7). Based on my quick skim, Hoskins et al. don't discuss how static stability can modulate vertical motions at all.

Thus, the extant literature and software out there seems about equally divided. I think we should leave it up to the user. If they call q_vector and don't supply a value for static_stability, that means they want the version of Q without division by $\sigma$, with units to match. If they call q_vector and do supply a value of static_stability, that means they want the version of Q that does divide by $\sigma$, with units to match. I haven't actually tested the code yet, but I think this is what currently happens. If this is correct, then my vote would be "not a bug".

I would agree with @dopplershift that just tacking on units to the default value of 1 is not a good solution, since a value of 1 with SI units would be a pretty implausible value for $\sigma$ in the real world.

Regarding this comment:

The error is in the static stability parameter that is the default value doesn't contain any units, even though the documentation states that the parameter should be a value with units, the default value violates that documentation!

The documentation I see says "static_stability (pint.Quantity, optional) – The static stability at the pressure level. Defaults to 1 if not given to calculate the Q-vector without factoring in static stability." This seems to match the current implementation as far as I can tell.