Using Virtual Temperature in CAPE/CIN calculations

[mgrover1]

Description Of Changes

Move toward using virtual temperature for the parcel profile instead of temperature

Closes #461
Closes #2093

Checklist

• Closes Add parcel path options #461
• Tests added
• Fully documented

[mgrover1]

While the tests are still failing (due to changing values) I would be interested to get feedback from @dopplershift / @jthielen on this approach... does this make sense? I tried to implement approaches discussed over in #461 ...

I would much rather get this fixed here instead of the downstream approach we have in ACT...

[dopplershift]

@mgrover1 I don't want to let this go too long without a response, but generally supportive of the approach. My first concern is making sure we get all the math right, which I haven't had a chance to sit down and make sure it all makes sense.

My other concern is making sure we can easily plumb this into all the places that need this knowledge for consistency.

[mgrover1]

@mgrover1 I don't want to let this go too long without a response, but generally supportive of the approach. My first concern is making sure we get all the math right, which I haven't had a chance to sit down and make sure it all makes sense.

My other concern is making sure we can easily plumb this into all the places that need this knowledge for consistency.

Thanks @dopplershift - I can put together some direct comparisons with the current method vs. the approach here, and look for a potential test case to make sure the math works out...

[mgrover1]

An Update

I put together a comparison of the implementation here compared to one currently in MetPy main

MetPy Main Branch

Surface-based CAPE: ~ 2329.4J/kg
Surface-based CIN: ~ -92.32 J/kg

This PR (using virtual temperature)

Surface-based CAPE: ~ 2720.4J/kg
Surface-based CIN: 0 J/kg

Is this right?

Graphically, this seems like the right approach, and takes into account the previous comments on this PR/related issues. What do people think? One thing to note here is that this will change the CAPE/CIN values in the testing suite, so we should either set the default value in the testing suite for whether to use the virtual temperature to be False, or the entire CAPE/CIN testing suite will need to be updated.

[Meteodan]

This still doesn't look correct. I'd have to take a look at the code, but from just looking at the plot above, the LCL should not be any different from the original uncorrected profile. It looks like you are using the virtual temperature at the surface in place of the regular temperature, and then lifting the parcel accordingly. This is leading to a higher LCL, and if so it's the wrong approach. (EDIT: yes, this appears to be what is happening - see my comments on the code below). Instead, the parcel should be lifted using the regular surface temperature as before, but then (and this is the crucial part), the resulting parcel temperature profile should be corrected using the virtual temperature. That is, compute the virtual temperature of the parcel given the parcel temp and dewpoint profile. then pass that profile into the CAPE calculation. Hopefully that makes sense.

Forgot the other part. When passing in the parcel virtual temp profile to the CAPE calculation, it should be compared against the environmental virtual temperature profile, not the environmental regular temperature profile. It looks like you are doing that though.

[Meteodan]

@mgrover1, I should be clear that I appreciate that you are taking the initiative here, and hope I'm not coming across as overly pedantic.

[dopplershift]

@Meteodan It's definitely appreciated to have the review of someone with intimate knowledge of what to expect.

[mgrover1]

@Meteodan agreed! Thank you for the helpful input here!!

[mgrover1]

@Meteodan I made those changes you suggested (I think) and this is the resultant SkewT diagram, highlighting the parcel profile and shaded CAPE... this doesn't look right either.

I think I am more confused than before...

and here is the code block changed

def _parcel_profile_helper(pressure, temperature, dewpoint,
                           assumptions=ParcelPathAssumptions()):
    """Help calculate parcel profiles.

    Returns the temperature and pressure, above, below, and including the LCL. The
    other calculation functions decide what to do with the pieces.

    """
    # Check that pressure does not increase.
    if not _check_pressure(pressure):
        msg = """
        Pressure increases between at least two points in your sounding.
        Using scipy.signal.medfilt may fix this."""
        raise InvalidSoundingError(msg)

    # Find the LCL
    press_lcl, temp_lcl = lcl(pressure[0], temperature, dewpoint)
    press_lcl = press_lcl.to(pressure.units)

    # Find the dry adiabatic profile, *including* the LCL. We need >= the LCL in case the
    # LCL is included in the levels. It's slightly redundant in that case, but simplifies
    # the logic for removing it later.
    press_lower = concatenate((pressure[pressure >= press_lcl], press_lcl))
    temp_lower = dry_lapse(press_lower, temperature)

    # Do the virtual temperature correction for the parcel below the LCL
    if assumptions.use_virtual_temperature:
        # Calculate the relative humidity, mixing ratio, and virtual temperature
        rh_lowest = relative_humidity_from_dewpoint(temperature, dewpoint)
        rv_lower = mixing_ratio_from_relative_humidity(press_lower[0], temperature, rh_lowest)
        temp_lower = virtual_temperature(temp_lower, rv_lower).to(temp_lower.units)
        temp_lcl = virtual_temperature(temp_lcl, rv_lower)

    # If the pressure profile doesn't make it to the lcl, we can stop here
    if _greater_or_close(np.nanmin(pressure), press_lcl):
        return (press_lower[:-1], press_lcl, units.Quantity(np.array([]), press_lower.units),
                temp_lower[:-1], temp_lcl, units.Quantity(np.array([]), temp_lower.units))

    # Establish profile above LCL
    press_upper = concatenate((press_lcl, pressure[pressure < press_lcl]))

    # Remove duplicate pressure values from remaining profile. Needed for solve_ivp in
    # moist_lapse. unique will return remaining values sorted ascending.
    unique, indices, counts = np.unique(press_upper.m, return_inverse=True, return_counts=True)
    unique = units.Quantity(unique, press_upper.units)
    if np.any(counts > 1):
        warnings.warn(f'Duplicate pressure(s) {unique[counts > 1]:~P} provided. '
                      'Output profile includes duplicate temperatures as a result.')

    # Find moist pseudo-adiabatic profile starting at the LCL, reversing above sorting
    temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units)
    temp_upper = temp_upper[::-1][indices]

    # Do the virtual temperature correction for the parcel above the LCL
    if assumptions.use_virtual_temperature:

        # Calculate the mixing ratio and virtual temperature
        rv_upper = mixing_ratio(saturation_vapor_pressure(temp_upper), press_upper)
        temp_upper = virtual_temperature(temp_upper, rv_upper).to(temp_lower.units)

    # Return profile pieces
    return (press_lower[:-1], press_lcl, press_upper[1:],
            temp_lower[:-1], temp_lcl, temp_upper[1:])

[Meteodan]

Don't correct the parcel profile until after the whole thing is computed. Here you are computing the lower part of the profile, then correcting it, then computing the upper part based on the corrected lower part, and then correcting the upper part. This is leading to the bizarre behavior near and above the LCL. Instead, compute both the lower and upper parts before the correction, and then apply the virtual temp correction to the whole thing.

[mgrover1]

Ahhh that makes more sense... I will give that a try! Sorry for the confusion here... thermo is slowly coming back :)

[sgdecker]

The classic paper on this subject contains a few soundings that could be used as additional test cases, assuming the Wyoming server has them.

[mgrover1]

The classic paper on this subject contains a few soundings that could be used as additional test cases, assuming the Wyoming server has them.

Thanks!! This is a great resource.

[Meteodan]

Think of it this way.... as far as the parcel is concerned, as it is lifted, its temperature will change based on the dry and moist adiabatic lapse rates. It's virtual temperature is just "along for the ride". The virtual temperature doesn't matter for the actual thermodynamics of the parcel behavior; it is only a diagnostic quantity that is a convenient way to incorporate the impact of water vapor on the density. The only time the virtual temperature matters here is when computing the CAPE because CAPE is the vertical integral of the buoyancy:

$B\approx -g\frac{\rho'}{\rho_{env}}\approx g\frac{T_v'}{T_{env}}=g\frac{T_{v,parcel}-T_{v,env}}{T_{v,env}}$

EDITED to correct formula for buoyancy (was in a hurry and forgot to multiply by g as well as the minus sign when using the density perturbation/environmental density)

[Meteodan]

This is why you should use the original parcel T and Td to compute its profile, including the dry adiabatic layer below the LCL, the LCL itself, and the moist adiabatic profile above. That's the true parcel temperature as it rises (assuming parcel theory is correct, of course). The virtual temperature of the parcel can then be computed at each level knowing the T and Td. It's this virtual temperature (ie. "corrected") profile that is passed (along with the environment's virtual temperature profile) to the CAPE routine to compute the CAPE with the "correction". FWIW, I don't like the word "correction" in this context. Really what we are doing is making a more accurate calculation of CAPE that includes the effect of water vapor on the density of the parcel and environment, and thus, the buoyancy.

[mgrover1]

Alrighty - @dopplershift this should be good for review. My only other question is whether we should update the cape_cin docstring to reflect that we are using virtual temperature

Tests are passing too 😄

[dopplershift]

Alrighty - @dopplershift this should be good for review. My only other question is whether we should update the cape_cin docstring to reflect that we are using virtual temperature

YES. Please do.

I'm going to go through the sounding in that paper with a fine-tooth comb to convince myself we're good now too. 😉

[mgrover1]

Alrighty - @dopplershift this should be good for review. My only other question is whether we should update the cape_cin docstring to reflect that we are using virtual temperature

YES. Please do.

I'm going to go through the sounding in that paper with a fine-tooth comb to convince myself we're good now too. 😉

@dopplershift here is a screenshot of the rendered docs updated in this PR...

[Meteodan]

Also, here is the more realistic looking Skew-T, showing the various temperature, virtual temperature, and dewpoint lines.

A student of mine pointed out something that I didn't notice at first about this plot. For the uncorrected parcel trace, the lapse rate below the LCL looks (dry) subadiabatic when it should be dry adiabatic....

[sgdecker]

For the uncorrected parcel trace, the lapse rate below the LCL looks (dry) subadiabatic when it should be dry adiabatic....

My guess is that the LCL point, while being calculated internally, is not included in the plot.

[mgrover1]

For the uncorrected parcel trace, the lapse rate below the LCL looks (dry) subadiabatic when it should be dry adiabatic....

My guess is that the LCL point, while being calculated internally, is not included in the plot.

Here is the plot with the LCL; I suspect it is more of an artifact of poor low-level resolution in the sounding. (1000, 1000, 950 hPa), with an LCL of 968 hPa.

The NASA Wallops sounding from the Doswell paper has an improved resolution, along with actual heights that we can check to ensure the parcel follows this idea of being dry adiabatic

One thing to note here is that MetPy does not assume a direct 9.8 degC/km for calculating the dry lapse rate below the LCL, but rather uses the potential temperature (see the function docs)

There is a nice example that digs into "The Mountain Problem", the classic atmospheric thermodynamics example, that explains this more in detail

• Link to mountain problem

From the AMS glossary,

So, the potential temperature at the surface is preserved up unto the LCL.

https://github.com/Unidata/MetPy/blob/main/src/metpy/calc/thermo.py#L1177

[Meteodan]

For the uncorrected parcel trace, the lapse rate below the LCL looks (dry) subadiabatic when it should be dry adiabatic....

My guess is that the LCL point, while being calculated internally, is not included in the plot.

Ok, that makes sense!

[Meteodan]

One thing to note here is that MetPy does not assume a direct 9.8 degC/km for calculating the dry lapse rate below the LCL, but rather uses the potential temperature (see the function docs)

Sure... that's pretty much the same thing, though. A parcel being lifted dry adiabatically will conserve its potential temperature, and its temperature will decrease at about 9.8 degC/km.

I think the culprit is what @sgdecker suggested above and what you showed by plotting the actual LCL: the plot doesn't have the necessary resolution so is interpolating between points above and below the LCL, leading to the apparent subadiabatic behavior.

[mgrover1]

@dopplershift - was there anything else here you think needs to be addressed?

[dopplershift]

I think we finally got there. Thanks for the hard work @mgrover1 !

And thanks to everyone for the detailed input and consideration!