Large negative CAPE values

[mkituseful]

What went wrong?

I have posted as a question at https://stackoverflow.com/questions/79556606/unrealistic-very-large-negative-cape-from-metpy-surface-based-cape-cin-and-mos, but thought it might need to be a bug report.

Operating System

Linux

Version

1.6.3

Python Version

3.12.8

Code to Reproduce

See the stackoverflow post.

Errors, Traceback, and Logs

[kgoebber]

Hi @mkituseful thanks for bringing this issue up. I think I have been able to track down what is going on with the calculation! In our cape_cin calculation we convert to using virtual temperature to reflect the best practice of computing those values. With the data you provided, the data go up to 1 hPa, which is way up in the thermosphere, typically well beyond where we test our calculations, especially for thermodynamic calculations that are best suited to the troposphere and lower stratosphere at most (typically >= 100 hPa). As it turns out, we run into some issues at the 2 hPa and 1 hPa levels when converting the parcel profile temperatures to virtual temperatures and that is the result of getting negative mixing ratios, which are non-physical. This is likely the result that those values are violating some assumptions when we are at so low of pressures.

import numpy as np
from metpy.units import units
import metpy.calc as mpcalc

pres = np.array([1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125, 100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1])
temp = np.array([27.8, 26.5, 24.2, 22.1, 20.0, 17.4, 15.2, 13.2, 11.4, 9.6, 7.8, 4.3, 0.5, -3.1, -7.0, -11.5, -16.2, -22.2, -30.3, -39.3, -49.3, -53.7, -57.7, -59.9, -59.7, -60.3, -61.7, -61.3, -57.8, -51.2, -48.2, -37.6, -32.5, -25.9, -16.3, -3.5, -4.9])
tdew = np.array([14.2, 13.2, 12.6, 11.9, 11.2, 10.7, 9.8, 9.1, 7.6, 5.1, 2.0, -2.3, -5.0, -10.6, -19.9, -23.5, -28.2, -32.8, -38.3, -48.7, -58.5, -62.2, -63.1, -68.3, -76.8, -83.7, -85.5, -87.2, -88.7, -91.3, -93.3, -96.6, -98.2, -99.7, -102.1, -103.8, -106.8])
pres = pres * units('hPa')
temp = temp * units('degC')
tdew = tdew * units('degC')

pressure_lcl, _ = mpcalc.lcl(pres[0], temp[0], tdew[0])
below_lcl = pres > pressure_lcl
parcel_mixing_ratio = np.where(below_lcl, mpcalc.saturation_mixing_ratio(pres, tdew),
                                mpcalc.saturation_mixing_ratio(pres, temp))

This produces the following output (scroll all the way to the right):

[0.01023027391291015 0.009825554242558883 0.009693778825702111 0.009504888997569188 0.00932425323700488 0.009276519130770934 0.00898787101602343 0.008832252120309662 0.010652107261577667 0.009737793395334728 0.00889693744260551 0.007467813543176856 0.0061236550707205045 0.00508594855990761 0.00412395726361982 0.003182723802102021 0.0024069812056159274 0.001618776159140922 0.0008828043976073874 0.00042285860332019493 0.00017139198690604418 0.00011393184534966621 7.86721195999433e-05 6.812321660454472e-05 8.153160443766854e-05 9.061851956935795e-05 9.453288358846891e-05 0.00014226016288342347 0.00031090460736564815 0.0011492673100951022 0.002435980145415396 0.015461054795193389 0.0378892791306624 0.10908176802041285 0.8362086411881933 -1.0786773875768485 -0.813214249559921]

If you limit the data to be more typical pressure levels (e.g., those greater than or equal to 100 hPa), then you get an expected outcome.

import numpy as np
from metpy.units import units
import metpy.calc as mpcalc

pres = np.array([1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 225, 200, 175, 150, 125, 100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1])
temp = np.array([27.8, 26.5, 24.2, 22.1, 20.0, 17.4, 15.2, 13.2, 11.4, 9.6, 7.8, 4.3, 0.5, -3.1, -7.0, -11.5, -16.2, -22.2, -30.3, -39.3, -49.3, -53.7, -57.7, -59.9, -59.7, -60.3, -61.7, -61.3, -57.8, -51.2, -48.2, -37.6, -32.5, -25.9, -16.3, -3.5, -4.9])
tdew = np.array([14.2, 13.2, 12.6, 11.9, 11.2, 10.7, 9.8, 9.1, 7.6, 5.1, 2.0, -2.3, -5.0, -10.6, -19.9, -23.5, -28.2, -32.8, -38.3, -48.7, -58.5, -62.2, -63.1, -68.3, -76.8, -83.7, -85.5, -87.2, -88.7, -91.3, -93.3, -96.6, -98.2, -99.7, -102.1, -103.8, -106.8])
pres = pres * units('hPa')
temp = temp * units('degC')
tdew = tdew * units('degC')

# CAPE/CIN
cape, cin = mpcalc.most_unstable_cape_cin(pres[pres >= 100*units.hPa],
                                          temp[pres >= 100*units.hPa],
                                          tdew[pres >= 100*units.hPa])
print(f"Most Unstable CAPE: {cape:.1f}, CIN: {cin:.1f}")
# Surface based
cape, cin = mpcalc.surface_based_cape_cin(pres[pres >= 100*units.hPa],
                                          temp[pres >= 100*units.hPa],
                                          tdew[pres >= 100*units.hPa])
print(f"Surface-based CAPE: {cape:.1f}, CIN: {cin:.1f}")

With the following output:

Most Unstable CAPE: 169.4 joule / kilogram, CIN: -65.6 joule / kilogram
Surface-based CAPE: 169.4 joule / kilogram, CIN: -65.6 joule / kilogram

Additionally, that comports with the quick graphic that I made to have a sanity check on the CAPE/CIN values

[DWesl]

Does #3758 fix your issue?

[mkituseful]

Hi, @kgoebber:

Thank you so much for spending time looking into this problem. According to your suggestions, I did further check on the more details. Initially I wondered how EL calculation wasn't able to locate the right level. Now, yes there were a couple levels near the top (< 5 mb) where mixing ratio are unphysical, which resulted in an intersection of the parcel temperature and environmental temperature, leading to a very high EL. With that fix, both unreasonably large negative and positive CAPE values disappeared.

Thanks a lot!