Analytic LCL

[jrleeman]

@brianmapes pointed out this paper: http://romps.org/pubs/pubs-2016-lcl.html that has an analytic LCL. The code on the author's webpage has some syntax/import issues. It would be interesting to compare this to our LCL implementation in terms of results and speed and possibly implement it in MetPy.

[dopplershift]

Here's a link to the same article, but to make more clear, it is a publication in an AMS journal: https://journals.ametsoc.org/doi/abs/10.1175/JAS-D-17-0102.1

[xebadir]

I've run this up against the Stull [2000] approximation and the iterative solver you guys use - its reasonably spritely and most definitely the most accurate approach (from the theoretical standpoint), but would probably require some playing under the hood to utilize existing metpy functions/structure. Part of the issue is it also handles LDL and minimum LCL height which adds to the complexity. Will post a comparison profile against current benchmarks.

[jthielen]

I've played around with this recently in my coursework, and I wrote up a new Python implementation for it using pint quantities. While I would still need to do some more robust tests, it looks like it is 2-3x faster than MetPy's current fixed-point solver formulation, and gives results within about 0.5 hPa and 0.05 K. Would there be interest in replacing MetPy's implementation with this? Or (at the risk of opening up a can of worms) having it as an alternative calculation?

[jthielen]

Also, in relation to #1199, it does work for grids.

[zbruick]

It would be great to at least see it! I think this is accurate enough for all practical purposes. If we gain performance too, all the better. Have you run it against the full test suite to catch edge cases?

[jthielen]

No, I haven't yet, but I'll plan on doing that soon. I'm also going to test the constants against those in MetPy, and try to see if any floating point issues arise given how many operations are used. Here's what I've been using so far:

from scipy.special import lambertw
from metpy.units import units
from metpy.calc import relative_humidity_from_dewpoint, specific_humidity_from_dewpoint

R_a = 287.04 * units('J/kg/K')
c_va = 719 * units('J/kg/K')
c_vv = 1418 * units('J/kg/K')
p_trip = 611.65 * units('Pa')
T_trip = 273.16 * units('K')
E_0v = 2.3740e6 * units('J/kg')
R_v = 461 * units('J/kg/K')
c_vl = 4119 * units('J/kg/K')
c_pa = c_va + R_a
c_pv = c_vv + R_v

def exact_lcl(pressure, temperature, dewpt):
    rh = relative_humidity_from_dewpoint(temperature, dewpt)
    q_v = specific_humidity_from_dewpoint(dewpt, pressure)
    R_m = (1-q_v) * R_a + q_v * R_v
    c_pm = (1 - q_v) * c_pa + q_v * c_pv
    a = c_pm / R_m + (c_vl - c_pv) / R_v
    b = -(E_0v - (c_vv - c_vl) * T_trip) / (R_v * temperature)
    c = b / a
    T_lcl = c / lambertw(rh**(1/a) * c * np.exp(c), k=-1).real * temperature
    p_lcl = pressure * (T_lcl / temperature)**(c_pm / R_m)
    return p_lcl, T_lcl

[zbruick]

I'm sure @dopplershift has thoughts on this too, but I think it would be a good replacement if we can make it happen and make sure we hit all of the edge cases we have in the test suite.

[dopplershift]

I don't know that I want to start the alternatives with this, but I think it would be reasonable to replace our existing with this. I actually like the idea that this came along after our original implementation, and then updated to use the latest science.

I do want to read through and try to understand the thermo here, but I have no qualms about moving the calculation.

[jthielen]

Here are the results in testing it against the thermo tests: https://gist.github.com/jthielen/9db772f0681bfee95b473a163669ca23

In short, there are many failed tests due to changes in results, and unfortunately the results are rather sensitive to the set of constants used (for instance, test_lcl gives 867.9171544080796 hPa with the Romps (2017) formulas and MetPy constants, and 863.9110172044186 hPa with the Romps (2017) formulas and constants, compared to the currently expected result of 864.761 hPa).

And so, I'd be hesitant in moving forward until there can be more discussion about what the "true" results should be.

[zbruick]

Thanks for providing the results. I agree that we need to nail down the science side of things before moving forward. There are also some edge cases that look like the analytic solution is returning a LFC now that were previously designated as NaNs, so I'm curious to see what's happening in those cases too.

[dopplershift]

@jthielen 🤢

Thanks for looking into it. There's definitely a lot more to be figured out before we can move on this.

[DLFinney]

FYI, David Romps is happy to help if someone wants to progress this... https://github.com/Unidata/MetPy/discussions/3029#discussioncomment-5776301

[robwarrenwx]

Hi MetPy team. First time contributor here. I'm a research scientist at the Australian Bureau of Meteorology. My team recently built a Python library for calculating a wide array of convective parameters (similar to SHARPpy but much faster) for use by our forecasters. We use the Romps (2017) non-iterative method for our LCL calculations, but with a slightly different formulation. Specifically, we use Cpv (the specific heat of water vapor at constant pressure) in place of Cvv (the specific heat of water vapor at constant volume) and Lv0 (the difference in specific enthalpy between water vapour and liquid water at the triple point) in place of Ev0 (the difference in specific internal energy between water vapour and liquid water at the triple point). The two sets of parameters are related as follows: Cpv = Cvv + Rv Lv0 = Ev0 + Rv * T0 where T0 is the triple point temperature (273.16 K). I've included the LCL function below for reference. We get almost perfect agreement between this formulation and an "exact" iterative calculation (using Newton's method). Note that we have a few checks in there to (a) handle (super)saturated initial conditions and (b) ensure that the LCL is not below (i.e. at a higher pressure than) the initial level, which can occur when the lifted parcel is close to saturation due to inaccuracies in the Lambert-W solver.

I suspect that part of the reason for the inconsistencies identified by @jthielen is that currently MetPy effectively ignores the variation of the specific gas constant R and the specific heat at constant pressure Cp with moisture. To get results consistent with the Romps formulation, you would need to replace the constant kappa in the _lcl_iter function with a locally defined "moist" kappa, given by kappa_m = kappa * (1 + w * (Rv / Rd)) / (1 + w * (Cpv / Cpd)) Note that if you're going to allow R and Cp (and therefore kappa) to vary with moisture here, you probably need to do so elsewhere to ensure that all calculations are thermodynamically consistent. For example, in the definition of potential temperature and the dry and moist lapse rates. For the moist lapse rate, you simply need to retain the term b from Eq. 8 of Bakhshaii and Stull (2013), which is identical to the quotient in the equation for kappa_m above.

It's also worth noting that Romps' LCL equations implicitly account for the variation of the latent heat of vaporization with temperature, which is also currently neglected in MetPy. This is a less valid assumption than constant R and Cp, and can have a big impact on thermodynamic quantities such as CAPE, as discussed, for example, by Yano and Ambaum (2017). This should probably be addressed in a separate issue.

Finally, I'll note that Romps also has derived an exact equation for the dewpoint temperature (Romps 2021), based on the same set of assumptions used for the LCL calculation (known as the Rankine-Kirchhoff approximations). It might be worth implementing this approach as a more precise alternative to the Bolton formulation. If you do this, you should also consider implementing the analytical formula for saturation vapor pressure used by Romps and derived separately by Ambaum (2020). Again, happy to raise this as a separate issue.

def lift_cond_level(p, T, q):
    """
    Computes pressure and parcel temperature at the LCL.

    Uses equations from Romps (2017).

    Romps, D.M., 2017. Exact expression for the lifting condensation
        level. Journal of the Atmospheric Sciences, 74, 3033-3057.

    Args:
        p: pressure (Pa)
        T: temperature (K)
        q: specific humidity (kg/kg)

    Returns:
        p_lcl: pressure at the LCL (Pa)
        T_lcl: temperature at the LCL (K)

    """

    # Compute relative humidity
    RH = rel_hum(p, T, q)
    if RH >= 1.:  # saturated conditions
        return p, T  # use initial values for LCL

    # Compute effective gas constant and specific heat
    Rm = eff_gas_const(q)
    Cpm = eff_spec_heat(q)

    # Set constants (Romps 2017, Eq. 22d-f)
    a = cpm / Rm + (Cpl - Cpv) / Rv
    b = -(Lv0 + (Cpl - Cpv) * T0) / (Rv * T)
    c = b / a

    # Compute temperature at the LCL (Romps 2017, Eq. 22a)
    fn = RH ** (1. / a) * c * np.exp(c)
    W = lambertw(fn, k=-1, tol=0.0001).real   # scipy function
    T_lcl = c * (1. / W) * T

    # Compute pressure at the LCL (Romps 2017, Eq. 22b)
    p_lcl = p * (T_lcl / T) ** (Cpm / Rm)

    # Check that LCL is not below initial level
    if p_lcl > p:  # inaccuracies in lambertw can cause this for q~=qs
        p_lcl = p
        T_lcl = T

    return p_lcl, T_lcl

[dopplershift]

@robwarrenwx @DLFinney Thanks the detailed analysis and bringing this up. I'm pretty convinced scientifically this is a robust, solid approach for us to move towards. The items I see:

1. How do any true constants used here compare to the values we now use. We tweaked some slightly with MetPy 1.0, and I'm really hoping we don't have to grapple with that again.
2. Updating tests is a going to be a pain, but doable. Need to make sure any tests for corner cases are still actually catching that corner (e.g. the missing LFC test mentioned above)
3. Curious, benchmark-wise, how this compares--especially important as I'm focused on performance. This might be a win-win here, but I want some data.

[dopplershift]

Romps (2021) for dewpoint is interesting (and we'd need to use that for saturation_vapor_pressure as well), but that is definitely best considered separately. Same for Romps (2015) for theta_e as well.

[robwarrenwx]

Hi @dopplershift - pleased to e-meet you!

Regarding the constants, below I've compared the relevant MetPy constants with those from Romps (2017). Note that I have assumed that Cvl ≈ Cpl ≡ Cl and Cvi ≈ Cpi ≡ Ci; however, this isn't strictly true, particularly for liquid water at higher temperatures (see Fig. 1 in Ambaum (2020)). Romps uses Cvl and Cvi whereas MetPy uses Cpl and Cpi. Furthermore, Romps optimized the values of Cvl and Cvi (along with Rv) to minimize the fractional difference between the saturation vapor pressures given in their Eqs. (4) and (5) and the saturation vapor pressures given by Wagner and Pruß (2002) and Wagner et al. (2011) based on laboratory data. This explains the larger percentage differences seen for Cl and Ci in the table. All the other constants are within ~1%. I'd say you're probably safe to stick with the constants you have. Ultimately, most of these quantities aren't truly constant so there are a range of values that could reasonably be use. This also means that a precision of 4 or 5 significant figures is plenty. For our Convective Parameters code, we use Rd = 287.04, Rv = 461.52, Cpd = 1005.7, Cpv = 1884.4, Cpl = 4219.9, and Lv0 = 2501000. The values of Cpv and Cpl are taken from Tabe 13.1 of Wagner and Pruß (2002) and are valid at the triple point (273.16 K). We don't currently consider ice processes so haven't specified Cpi or Lf0.

Constant	MetPy	Romps	Diff. (%)
Rd	287.047...	287.04	-0.00261
Rv	461.523...	461	-0.113
Cpd	1004.666...	1006.04	0.137
Cvd	717.618...	719	0.192
Cpv	1860.078...	1879	1.02
Cvv	1398.554...	1418	1.39
Cl	4219.4	4119	-2.38
Ci	2090	1861	11.0
Lv0	2500840	2499927	-0.0365
Lf0	333700	333700	0

[robwarrenwx]

@dopplershift

As far as efficiency goes, I just did a quick test comparing my implementation of Romps (2017) with the MetPy LCL function for temperatures between -50 and 50 degC and relative humidites between 1 and 99 % at a pressure of 1000hPa. The Romps function is generally around 20-30 times faster. More importantly, since the SciPy Lambert-W function can handle arrays, the Romps implementation could easily be vectorized, allowing for the efficient calculation of LCL pressure and temperature over multidimensional arrays.

With the above-noted tests, I found differences in LCL pressure and temperature of up to 2 hPa and 0.15 degC, respectively. These don't reflect errors in the Romps method, as his equations are exact (although small errors may be introduced in the Lambert-W solver). Nor do they reflect different constants, as I used the values from MetPy. Instead, I believe they are caused by other approximations in MetPy, such as the assumption of constant R and Cp and the use of Bolton's formula for saturation vapour pressure (rather than the analytical expression used by Romps and presented separately by Ambaum 2020).

To ensure that moist thermodynamics is treated consistently, I would recommend the following changes:

1. Implement functions to calculate Rm, the effective gas constant for moist air (Romps 2017, Eq. 16) and Cpm, the effective specific heat at constant pressure for moist air (Romps 2017, Eq. 17).
2. Implement functions to calculate the latent heats of vaporization and freezing as a function of temperature (Ambaum 2020, Eq. 15 and 18).
3. Update the saturation vapour pressure function to use the analytical equation from Ambaum (2020) (his Eq. 13; equivalent to Eq. 4 from Romps 2017). Consider implementing an additional function to calculate the saturation vapour pressure over ice (Ambaum 2020, Eq. 17; equivalent to Eq. 5 from Romps 2017).
4. Update the LCL function to use the analytical equations (Romps 2017, Eq. 22). Consider also adding a function to calculate the lifted deposition level (LDL; Romps 2017, Eq. 23). Note that the equations for the height of the LCL and LDL (Eq. 22c and 23c) should not be used as they are inexact. It is better to get the height via interpolation from the LCL/LDL pressure.
5. Implement the analytical function for dewpoint temperature from Romps (2021) (his Eq. 5 and 6). Consider also adding a function to calculate the frost-point temperature (Romps 2021, Eq. 7 and 8).

I propose that we raise a separate issue to address 1-3. With these changes, modifications will be needed to other functions that currently ignore the moisture dependence of R and Cp (such as potential_temperature) and/or the temperature dependence of Lv (such as moist_lapse). New unit tests may also be needed. The LCL function should only be updated once these changes are implemented. As noted above, use of the Romps (2021) functions for dewpoint (and, optionally, frost-point) temperature should also be raised as a separate issue.