MetPy计算ERA5数据出现超大负CAPE问题排查求助

使用ERA5再分析气压层数据的单个格点，调用MetPy的surface_based_cape_cin()和most_unstable_cape_cin()函数时，返回了极大的负CAPE值。手动计算能得到合理的CAPE结果，初步排查发现MetPy的cape_cin函数可能因未识别出EL（平衡高度），而将积分范围延伸至数据顶层（1hPa），最终积分得到超大负值，但不确定是否与find_intersections()函数有关。

示例代码及输出 #

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, temp, tdew)
print(f"Most Unstable CAPE: {cape:.1f}, CIN: {cin:.1f}")
# Surface based
cape, cin = mpcalc.surface_based_cape_cin(pres, temp, tdew)
print(f"Surface-based CAPE: {cape:.1f}, CIN: {cin:.1f}")

输出结果：

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

手动计算代码及结果 #

# Parcel profile 
prof = mpcalc.parcel_profile(pres, temp[0], tdew[0]).to('degC')
print("profile:", "[" + ", ".join(f"{x:.1f}" for x in prof.magnitude) + "]", prof.units)

# Find LFC
lfc_pres, lfc_temp = mpcalc.lfc(pres, temp, tdew, parcel_temperature_profile=prof, which='bottom')
print(f"LFC Pressure: {lfc_pres:.1f}, LFC Temperature: {lfc_temp:.1f}")

# Find EL
el_pres, el_temp = mpcalc.el(pres, temp, tdew, parcel_temperature_profile=prof, which='top')
print(f" EL Pressure: {el_pres:.1f},  EL Temperature: {el_temp:.1f}")

# CAPE calculation between LFC and EL
if lfc_pres is not np.nan and el_pres is not np.nan:
    # Find indices for LFC and EL
    lfc_idx = np.where(pres >= lfc_pres)[0][-1]
    el_idx = np.where(pres <= el_pres)[0][0]

    # Integrate CAPE (positive area only)
    cape = 0.0
    for i in range(lfc_idx, el_idx):
        if prof[i] > temp[i]:
            delta_ln_p = np.log(pres[i+1].magnitude / pres[i].magnitude)
            T_env = temp[i].to('K').magnitude        # In Kelvin
            T_diff = prof[i].magnitude - temp[i].magnitude
            cape += 9.8 * (T_diff / T_env) * delta_ln_p
    cape = (-cape * 287 * T_env) * units('J/kg')
else:
    cape = 0 * units('J/kg')

# CIN from MetPy
cin = mpcalc.cape_cin(pres, temp, tdew, prof)[1]

print(f"custom CAPE: {cape:.1f}, CIN: {cin:.1f}")

输出结果：

profile: [27.8, 25.6, 23.4, 21.2, 18.9, 16.5, 14.1, 11.7, 10.2, 8.9, 7.6, 4.8, 1.6, -2.0, -6.0, -10.6, -16.1, -22.5, -30.3, -39.7, -50.9, -57.3, -64.4, -72.2, -80.8, -90.5, -101.8, -118.4, -132.6, -151.7, -165.0, -184.4, -193.0, -200.4, -210.2, -217.1, -227.2] degree_Celsius
LFC Pressure: 734.9 hectopascal, LFC Temperature: 6.8 degree_Celsius
EL Pressure: 349.0 hectopascal, EL Temperature: -30.5 degree_Celsius
custom CAPE: 1072.7 joule / kilogram, CIN: -66.3 joule / kilogram

问题分析与解决 #

• 从手动计算结果可知，该案例中LFC和EL均存在，环境廓线与气块廓线之间有明确的正浮力区域，对应合理的CAPE值（1072.7 J/kg）。
• MetPy内置函数返回负CAPE的核心原因：函数默认在未检测到EL时，会将积分范围延伸至数据顶层。而本案例中，气块廓线在EL（349hPa）以上远低于环境温度，积分该区域会产生巨大的负贡献，完全抵消了正CAPE的数值。
• 可行的解决方式：
  1. 手动计算CAPE：先通过mpcalc.lfc()和mpcalc.el()确定积分范围，只对LFC到EL之间的正浮力区域积分。
  2. 截取气压层：调用MetPy函数前，先截取LFC到EL对应的气压层数据，再计算CAPE。
  3. 检查函数参数：确认MetPy是否提供控制积分范围的参数，避免积分到数据顶层。

内容的提问来源于stack exchange，提问作者Xingang Fan