Use parcel prof for parcel mix ratio

[ahijevyc]

This PR has 2 bug fixes for deriving parcel mixing ratio in metpy.calc.thermo.cape_cin.

1. Below the LCL, set the parcel mixing ratio to a constant value equal to the mixing ratio at the surface. Do not use the pressure and dewpoint above the surface.
2. Above the LCL, use the parcel temperature to derive parcel mixing ratio, not the environmental temperature.

I noted this because this led to a negative mixing ratio when analyziing a test case sounding for CM1.

input_sounding_jordan_allmean: West Indies annual mean (Jordan, 1958, J. Meteor., p. 91)

in which the top level of the sounding is 40000 m, the potential temperature is 1095 K, the pressure is 10 hPa, and the temperature is 295 K. Plugging in 10 hPa and 295 K to metpy.calc.thermo.mixing_ratio, one gets a mixing ratio of -1.006 kg/kg. A negative mixing ratio is unphysical. Now I'm not sure if the mixing_ratio function has a bug also, but the problem goes away when you replace the environment temperature temperature with the much lower parcel temperature parcel_profile here link to code.

Description Of Changes

1. Change pressure, dewpoint to pressure[0], dewpoint[0] below LCL.
2. Change temperature to parcel_profile above LCL.
3. Update expected CAPE/CIN in tests and function docstrings.
4. Remove test for "sensitive" sounding calc.lfc( ) bug report. #902, which no longer exhibits positive CAPE.

[CLAassistant]

All committers have signed the CLA.

[DWesl]

Calculating the parcel mixing ratio on a pseudo-adiabat should definitely be done using the parcel temperature rather than the environment temperature.

mpcalc.mixing_ratio takes partial pressure and total pressure in that order, not total pressure and temperature. It could easily produce a negative mixing ratio if one switched the partial and total pressures, as then the partial pressure of water vapor would be larger than the total pressure, requiring a negative partial pressure of dry air.

mpcalc.saturation_mixing_ratio takes total pressure and temperature, and is what is called in the changes you made. It could also easily produce a negative mixing ratio if the total pressure is below the saturation vapor pressure for the given temperature, which I suspect might happen in the upper stratosphere. I think attempting to replicate this physically would produce a blob of pure water vapor, with mixing ratio infinity. The documentation should probably note the possibility and meaning, though that might be out-of-scope for this PR.

[ahijevyc]

@DWesl - Thanks for the quick response about the parcel temperature bug, and the thoughts about negative mixing ratios.

I meant to say I plugged in 10hPa and 295K to saturation_mixing_ratio, not mixing_ratio.

mpcalc.saturation_mixing_ratio(10*units.hPa, 295*units.K)
-1.0061825918566147 dimensionless

Sorry about the confusion.

I agree that addressing negative mixing ratio could be handled in another PR, and is beyond the scope of this PR. Perhaps the mixing_ratio function should return NaN and/or warn the user when the partial pressure exceeds the total pressure, instead of returning a negative mixing ratio.

[dopplershift]

Yeah, it seems like there's a legitimate bug fix here in regards to the virtual temperature correction. It looks like the changes you made for some reason include the changes in #3735; can you rebase onto main? I'm also guessing there are some tests that will need (likely minor) updates to pass with this change.

On the mixing ratio, yes the problem with the example in question is that the saturation vapor pressure is 26 hPa at 295 K, which is quite a bit larger than the given total pressure of 10 hPa; essentially at this point you can't actually condense because liquid water is boiling. I suspect the correct behavior in this case would be to return NaN (though I'd like to check some other tools), but that should be a separate PR.

[dopplershift]

There's still 1 extra commit (b55773f) on there.

[dopplershift]

I can do the rebase and force push it to your branch if you like. Then just need test updates before merging in.

[ahijevyc]

Thanks for the offer! I just did it. Think I just needed to eat some lunch first.
The tests will require some changes because the cape will change.

[ahijevyc]

Could I make the changes to the expected cape in the test files and work them into the PR?

[dopplershift]

Yes, please make the test changes in this PR so that all the tests pass before we merge the changes.

[ahijevyc]

I fixed the expected cape and cin in all the tests except test_sensitive_sounding. This one may require adjustments to the sounding temperature profile because it was supposed have some positive area, but with the new code, there is none.

https://github.com/Unidata/MetPy/blob/489b87cc707a1974b1579cf050f61f931e4fac7a/tests/calc/test_thermo.py#L746C1-L766C1

def test_sensitive_sounding():
    """Test quantities for a sensitive sounding (#902)."""
    # This sounding has a very small positive area in the low level. It's only captured
    # properly if the parcel profile includes the LCL, otherwise it breaks LFC and CAPE
    p = units.Quantity([1004., 1000., 943., 928., 925., 850., 839., 749., 700., 699.,
                        603., 500., 404., 400., 363., 306., 300., 250., 213., 200.,
                        176., 150.], 'hectopascal')
    t = units.Quantity([25.1, 24.5, 20.2, 21.6, 21.4, 20.4, 20.2, 14.4, 13.2, 13., 6.8, -3.3,
                        -13.1, -13.7, -17.9, -25.5, -26.9, -37.9, -46.7, -48.7, -52.1, -58.9],
                       'degC')
    td = units.Quantity([21.9, 22.1, 19.2, 20.5, 20.4, 18.4, 17.4, 8.4, -2.8, -3.0, -15.2,
                         -20.3, -29.1, -27.7, -24.9, -39.5, -41.9, -51.9, -60.7, -62.7, -65.1,
                         -71.9], 'degC')
    lfc_pressure, lfc_temp = lfc(p, t, td)
    assert_almost_equal(lfc_pressure, 952.8445 * units.mbar, 2)
    assert_almost_equal(lfc_temp, 20.94469 * units.degC, 2)


    pos, neg = surface_based_cape_cin(p, t, td)
    assert_almost_equal(pos, 0.106791 * units('J/kg'), 3)
    assert_almost_equal(neg, -282.620677 * units('J/kg'), 3)

This test was designed around issue #902.

[ahijevyc]

I was so focused on fixing the parcel mixing ratio high in the atmosphere that I forgot to consider below the LCL. The code also needs to use the parcel profile below the LCL, not the environmental profile. So I changed the logic below the LCL to use pressure[0] and dewpoint[0] to derive the parcel mixing ratio below the LCL (not environment pressure and dew point).

    # The mixing ratio of the parcel comes from the dewpoint below the LCL, is saturated
    # based on the temperature above the LCL
    parcel_mixing_ratio = np.where(below_lcl, saturation_mixing_ratio(pressure[0], dewpoint[0]),
                                   saturation_mixing_ratio(pressure, parcel_profile))

[dopplershift]

That seems correct as well. I'll see if I can find some values for that one sounding.

[DWesl]

In the spirit of the other changes, should the LCL be derived using parcel_profile[0] rather than temperature[0]? Frequently they're the same, but MLCAPE and MUCAPE would use other values.

[dopplershift]

@DWesl I'm not convinced that's correct. The use of virtual temperature in CAPE is to get a complete look at the buoyancy of a rising air parcel relative to the environment; virtual temperature accounts for the contribution of water vapor to that buoyancy. There is no buoyancy process in the LCL calculation, but rather you're lifting a parcel, assuming adiabatic dry processes with conserved contained water vapor (constant mixing ratio). That is, the LCL denotes the height at which condensation (and a shift to moist pseudo adiabatic processes) occurs--overall buoyancy plays no in determining that height, only the starting temperature and moisture content.

[DWesl]

parcel_profile is not converted from temperature to virtual temperature for another three lines.

I suspect you have more experience with such things than I do, so if you say Mixed-Layer CAPE and Most-Unstable CAPE use the environment temperature and mixing ratio at the surface to find where the mixed-layer or most-unstable parcel would form clouds, or that this function is entirely unsuited for MLCAPE/MUCAPE, I will concede to your experience. (With #3572 not in yet, MLCAPE will be tricky, and MUCAPE may well be covered by applying the current process to soundings starting on successive levels above the surface and taking the largest).

[ahijevyc]

Just to clear things up, I think @DWesl was suggesting that we use parcel temperature, as opposed to environment temperature, to calculate LCL (nothing to do with virtual temperature).

[dopplershift]

@ahijevyc Thanks for that, I got completely confused. My bad. 🐑

[dcamron]

I think using parcel_profile is more correct in support of arbitrary parcel profiles. In practice, when using profiles calculated by metpy.calc.parcel_profile, it won't make a difference. most_unstable_cape_cin and mixed_layer_cape_cin already massage input environment and parcel profiles such that this won't make a difference either. Either way won't account for assumptions made in parcel profiles calculated outside of metpy.calc.

However, this reveals to me that our parcel LCL > parcel mixing ratio > parcel virtual temperature and our LFC calculations assume that the parcel and the environment have the same initial moisture properties, even though cape_cin supports arbitrary parcel temperature profiles, which may have been calculated with differing moisture properties.

I don't think that's something we're going to fix cleanly in time for this release, so I will open a new issue. In the mean time, feel free to go ahead and update the LCL input to be the parcel profile temperature. Stick around for future discussion 😅

[ahijevyc]

The PR still fails one of the tests. It fails the "sensitive sounding", which was finely tuned to have a small but non-zero CAPE when the LCL was included in the profile and zero CAPE when the LCL was not included. With the PR change the CAPE is zero when the LCL is included.

I'm not sure this test is needed since we always include the LCL.

[ahijevyc]

With the bug fix implemented in this PR, one of the test soundings that had a small amount of positive CAPE ended up having zero CAPE. It was not practical to recreate a test sounding for which positive CAPE depended on whether the LCL was included or not.

Therefore, I removed the thermodynamic test under tests/calc/test_thermo/test_sensitive_sounding to allow the PR to continue the CI process.

After updating CAPE and CIN in a few docstring examples, the only failing tests now are in Build Docs / Linux X.y. Copilot suggests this may be an intermittent problem with the matplotlib web site.

[dopplershift]

Thanks for the fix! I agree, removing the test at this point seems best.