Problem with metpy.calc.lifted_index

[fuliii1]

Hi Team,

I test metpy.calc.lifted_index calculation with University of Wyoming Upper air observations data. In this database the lifted index is provided but if I calculate lifted index with Metpy I get different result and sometimes runtime warnings.

In the code example I used this database: https://weather.uwyo.edu/cgi-bin/sounding?region=europe&TYPE=TEXT%3ALIST&YEAR=2023&MONTH=07&FROM=1412&TO=1412&STNM=12843

from metpy.calc import dewpoint_from_relative_humidity, lifted_index, parcel_profile
from metpy.units import units

p = [1002.0, 1000.0, 993.0, 925.0, 850.0, 846.0, 843.0, 826.0, 823.0, 766.0,
    748.0, 734.0, 723.0, 711.0, 700.0, 686.0, 668.0, 658.0, 637.0, 632.0,
    579.0, 568.0, 550.0, 549.0, 511.0, 500.0, 479.0, 439.0, 420.0, 400.0,
    373.0, 356.0, 300.0, 284.0, 280.0, 273.0, 267.0, 250.0, 242.0, 241.0,
    239.0, 230.0, 215.0, 213.0, 212.0, 207.0, 204.0, 200.0, 196.0, 195.0,
    163.0, 150.0, 140.0, 135.0, 131.0, 127.0, 115.0, 112.0, 108.0, 100.0,
    93.0, 91.0, 87.1, 87.0, 81.0, 79.0, 77.0, 73.0, 72.7, 70.0, 69.0,
    67.0, 66.0, 64.0, 62.0, 54.6, 52.0, 51.0, 50.0, 47.0, 44.0, 42.0,
    41.0, 38.0, 37.0, 36.4, 36.0, 32.0, 31.0, 30.0, 29.0, 28.3, 28.0,
    26.0, 25.0, 23.0, 22.8, 20.0, 19.0, 18.0, 17.0, 16.4, 15.6, 15.0,
    14.0, 11.0, 10.0, 9.2, 9.0, 8.6
] * units.hPa

T = [28.2, 27.0, 25.4, 19.4, 12.8, 12.3, 12.0, 10.8, 11.2, 6.4, 
     4.8, 3.6, 4.2, 4.0, 3.8, 4.0, 2.2, 1.9, 1.2, 0.8, -3.5, -4.3,
     -5.7, -5.8, -9.9, -10.9, -12.7, -16.9, -18.3, -21.1, -25.2,
     -27.9, -38.5, -41.7, -42.6, -44.3, -45.5, -49.5, -51.7, -51.7,
     -52.3, -54.3, -57.5, -57.9, -58.1, -57.1, -56.5, -57.5, -58.5,
     -58.3, -52.3, -52.9, -53.9, -54.5, -54.9, -55.4, -56.9, -57.3,
     -56.4, -54.5, -57.0, -57.8, -59.3, -59.3, -59.0, -58.9, -58.8,
     -58.5, -58.5, -57.5, -57.2, -56.6, -56.3, -55.6, -55.0, -52.3,
     -52.9, -53.1, -53.3, -53.3, -53.4, -53.4, -53.4, -53.5, -53.5,
     -53.5, -53.3, -51.4, -50.8, -50.3, -49.4, -48.7, -48.8, -49.2,
     -49.4, -49.9, -49.9, -47.9, -47.3, -46.6, -45.9, -45.5, -45.9,
     -45.0, -43.4, -37.7, -35.5, -34.5, -34.6, -34.7
] * units.degC

rh = [.42, .45, .44, .60, .84, .86, .88, .86, .67, .76,
      .79, .82, .41, .23, .13, .6, .23, .17, .8, .27,
      .15, .33, .11, .11, .37, .28, .10, .46, .37, .37,
      .29, .25, .24, .20, .23, .28, .58, .67, .73, .73,
      .48, .42, .56, .49, .46, .26, .18, .20, .20, .19,
      .2, .1, .1, .1, .1, .1, .1, .1, .1, .1, .1, .1, 
      .1, .1, .1, .1, .1, .1, .1, .1, .1, .1, .1, .1,
      .1, .1, .1, .1, .1, .1, .1, .1, .1, .1, .1, .1,
      .1, .1, .1, .1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0
] * units.dimensionless

Td = dewpoint_from_relative_humidity(T, rh)
prof = parcel_profile(p, T[0], Td[0])
lift_index = lifted_index(p, T, prof)
print(f'   Lifted Index: {lift_index:.2f}')

lifted index provided by University of Wyoming: Lifted index: 2.02 LIFT computed using virtual temperature: 1.72

But the calculated with Metpy is Lifted Index: [-0.25] delta_degree_Celsius

The runtime warning logs: C:\Users\fulopg\AppData\Roaming\Python\Python312\site-packages\metpy\calc\thermo.py:1396: RuntimeWarning: divide by zero encountered in log val = np.log(vapor_pressure / mpconsts.nounit.sat_pressure_0c) C:\Users\fulopg\AppData\Roaming\Python\Python312\site-packages\metpy\calc\thermo.py:1397: RuntimeWarning: invalid value encountered in divide return mpconsts.nounit.zero_degc + 243.5 * val / (17.67 - val) Lifted Index: [-0.25] delta_degree_Celsius

Please help to judge which is the accurate lifted index.

Thanks for you help, Gabor Fulop

[fuliii1]

Hi Team,

I've checked the lifted index calculation formula in your documentation

and a reference documentation https://unidata.github.io/MetPy/latest/api/references.html#doswellschultz2006

For me it seems there are different formula in the documentation and in the referenced scientific article.

Could you help me which are the correcnt formula.

Thanks, Gabor Fulop

[dopplershift]

We have the correct formula in MetPy's docs and implementation of lifted_index. I just verified this with the AMS Glossary and Wikipedia (and in fact Wyoming itself uses that definition). The original paper (Galway 1956) doesn't give a formula but states:

The Lifted Index [2], as used in the SELS Center, is defined as the temperature difference between the observed 500 mb temperature and the assumed 500 mb temperature of a mean parcel lifted from the modified lower 3000 foot layer next to the ground. Indices are negative for parcel temperatures which are warmer than the environment.

This text description matches what we have implemented in MetPy. As far as I can tell, that's just an error in the Doswell and Schultz (2006) paper.

As to the discrepancy, there are a couple things. First, if we're trying to match Wyoming's calculated values, I would use the dewpoint values they provide rather than using MetPy to calculate since they calculate dewpoint differently. We can use Siphon to grab that data easily and reproduce your results above:

from datetime import datetime, UTC

import metpy.calc as mpcalc
from metpy.units import pandas_dataframe_to_unit_arrays, units
from siphon.simplewebservice.wyoming import WyomingUpperAir

data = WyomingUpperAir.request_data(datetime(2023, 7, 14, 12, tzinfo=UTC), '12843')
data = pandas_dataframe_to_unit_arrays(data)

prof = mpcalc.parcel_profile(data['pressure'], data['temperature'][0], data['dewpoint'][0])
lift_index = mpcalc.lifted_index(data['pressure'], data['temperature'], prof)

This gives -0.3488455533987125 delta_degree_Celsius, similar to what you got but slightly different because this uses the provided dewpoint values.

The more important change is to note what Galway (1956) says above about what parcel to use, "mean parcel lifted from the modified lower 3000 foot layer next to the ground"; Wyoming's documentation indicates they use a mixed parcel from the surface to 500m, so let's do the same using MetPy:

# Calculate 500m mixed parcel
parcel_p, parcel_t, parcel_td = mpcalc.mixed_parcel(data['pressure'], data['temperature'],
                                                    data['dewpoint'], height=data['height'])

# Replace sounding temp/pressure in lowest 500m with mixed values
above = data['height'] > 500 * units.m
pres = concatenate([parcel_p, data['pressure'][above]])
temp = concatenate([parcel_t, data['temperature'][above]])

# Calculate parcel profile along all desired pressure values from our mixed starting values
mixed_prof = mpcalc.parcel_profile(pres, parcel_t, parcel_td)

# Calculate lifted index using our mixed profile
lift_index = mpcalc.lifted_index(pres, temp, mixed_prof)

This gives a value of 2.125725101166301 delta_degree_Celsius, which is more in line with what you see on Wyoming's site. Any remaining difference is likely due to differences in how the parcel profile is calculated (esp. differences in calculating moist pseudoadiabatic ascent).

[dopplershift]

I'm thinking we should probably document the common use of a mixed parcel within the lifted_index documentation.

[dopplershift]

Could also update the example in the docstring to use a mixed parcel (and eliminate the need for dewpoint_from_relative_humidity(T, rh) while we're at it).

[fuliii1]

Hi @dopplershift,

Thanks for the clarification.

[dopplershift]

Re-opening to capture need for docs.

[sgdecker]

We have the correct formula in MetPy's docs and implementation of lifted_index. I just verified this with the AMS Glossary and Wikipedia (and in fact Wyoming itself uses that definition). The original paper (Galway 1956) doesn't give a formula but states:

The Lifted Index [2], as used in the SELS Center, is defined as the temperature difference between the observed 500 mb temperature and the assumed 500 mb temperature of a mean parcel lifted from the modified lower 3000 foot layer next to the ground. Indices are negative for parcel temperatures which are warmer than the environment.

This text description matches what we have implemented in MetPy. As far as I can tell, that's just an error in the Doswell and Schultz (2006) paper.

But look at the next paragraph in Galway (1956) [my italics for emphasis]:

The Lifted Index is determined in the following manner. First, the mean mixing ratio observed in the 3000 feet next to the ground is determined graphically (equi-areal method). Then, where significant heating (or cooling) is expected, the actual sounding is modified through the lower 3000 feet. This modification is made by assuming a dry adiabatic temperature lapse rate (lower 3000 feet only) through the predicted afternoon maximum temperature.

So Doswell and Schultz aren't really wrong. In practice, there are indeed much more common parcel definitions used to compute lifted index (like various mean parcels, the most unstable parcel, etc.), although I note that nsharp does have a "forecast surface" parcel option. I am not advocating using the "forecast surface" parcel approach in any documentation, or even coming up with a "forecast surface" function, since nowadays we have things like gridded model forecasts that did not exist in 1956, and so you can calculate the LI in the spirit of Galway by taking a mean parcel from a model forecast sounding rather than using the ad-hoc "forecast surface" approach. My main point is just to acknowledge that the original Galway method isn't quite as simple as it may seem, and if for some reason you are restricting yourself to analyzing 12Z observed soundings, it is probably the best way.

Edit: Just noticed the order of the subtraction is wrong in Doswell and Schultz. There's no disputing that! I was focused on the subscripts.

[dopplershift]

Edit: Just noticed the order of the subtraction is wrong in Doswell and Schultz. There's no disputing that! I was focused on the subscripts.

Yeah, my problem wasn't the subscripts, it's the order. Especially when Galway says:

Indices are negative for parcel temperatures which are warmer than the environment. Thus the Lifted Index is similar to the Showalter stability index except for the determination of the level from which the parcel is lifted and the fact that the Lifted Index is a forecast index whereas the Showalter Index is an observed static index.

and Doswell and Schultz (2006) have the signs opposite from each other. 🤷‍♂️

[fuliii1]

Hi @dopplershift, @sgdecker,

Thanks for the investigation and clarification.