Most unstable CAPE and surface-based CAPE comparison in NCL and Python

[lsun]

Hello,

I am testing NCL functions (wrf_cape_2d, wrf_cape_3d; wrf_user_getvar) for the CAPE calculation. What I did is to take a WRF output file provided by NCL website and use these functions to calculate the most unstable CAPE (wrf_cape_2d or wrf_user_getvar) and surface based CAPE (wrf_cape_3d or wrf_user_getvar). It came to my attention that the value of most unstable CAPE is smaller than surface-based CAPE (see my attached figures). Also surprising is that the value of most unstable CAPE using wrf_cape_2d function is a little larger than the most unstable cape using wrf_user_getvar(a,"cape_2d"). [By contrast, the CAPE calculated by wrf_cape_3d is identical to the CAPE calculated by wrf_user_getva(a,"cape_3d")]. I wonder if any one can help share some insight why this happens. In particular, shouldn't we expect the most unstable CAPE to be equal or larger than surface-based CAPE?

Note that I also tested this using WRF_Python package and found that wrf.cape_2d (identical to NCL wrf_cape_2d) is also smaller than wrf.cape_3d surface-layer value ( identical to NCL wrf_cape_3d). Thus, I suspect that this might be a common issue for both NCL and Python.

Any comments and/or suggestions will be appreciated!

--
Lantao sun
Research Scientist
Department of Atmospheric Science
Colorado State University
Fort Collins, Colorado, USA
Lantao.Sun@colostate.edu

P.S.: Some details below regarding how I did the test:
1. Most unstable CAPE based on wrf_cape_2d: I simply follow the code for example 1 at NCL website:
begin
a = addfile("wrfout_d01_2005-08-28_00:00:00","r")

T = a->T
P = a->P
PB = a->PB
QV = a->QVAPOR
PH = a->PH
PHB = a->PHB
HGT = a->HGT
PSFC = a->PSFC

T = T + 300.
P = P + PB
tk = wrf_tk( P , T )
PH = PH + PHB
z = wrf_user_unstagger(PH,PH@stagger)
z = z/9.81

cinfo = wrf_cape_2d( P, tk, QV, z, HGT, PSFC, True )

mcape = cinfo(0,:,

; ....................plot scripts follow ....................
end

2. Most unstable CAPE based on wrf_user_getvar: I follow the code for example 2 at NCL website:
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
a = addfile("wrfout_d01_2005-08-28_00:00:00","r")
cape2d = wrf_user_getvar(a,"cape_2d",-1)
mcape = cape2d(0,:,

3. Surface-based CAPE based on wrf_cape_3d: similar to 1, but change wrf_cape_2d to wrf_cape_3d, and only keep the bottom layer for surface-based CAPE.
4. Surface-based CAPE based on wrf_user_getvar: similar to 2, but use "wrf_user_getvar("cape3d",-1)" and only keep the bottom layer for surface-based CAPE.

 

[epotter1]

Hi,

I realise this is a very old post, but have found the same problem. Did you ever solve it? The difference between wrf_cape_2d and wrf_user_getvar(a,"cape_2d") is particularly weird.

Thanks for any help!

 

[lsun]

Unfortunately, no. I would appreciate it if someone familiar with these functions could provide some help.

 

[epotter1]

I'm using wrf-python rather than NCL, but my problem was with slight differences in the input values, I think. I ended up looking at the g_cape.py module and comparing all the input variables to find the difference.

Sorry for the horrible formatting and the python rather than NCL code, but in case it's helpful, the two seem to match up using the following:

<begin code>

from netCDF4 import Dataset
import wrf
import xarray as xr

file=Dataset('wrfout_file')
mytimeidx=3
test=xr.open_dataset(file)

### way 1
cape_2d_way1=wrf.getvar(file,'cape_2d',timeidx=mytimeidx)

### way 2
pres_hpa=(test.P[mytimeidx,:,:,:]+test.PB[mytimeidx,:,:,:])/100
qv=test.QVAPOR[mytimeidx,:,:,:]
height=wrf.getvar(file,'z',timeidx=mytimeidx)
terrain=test.HGT[mytimeidx,:,:]
psfc_hpa=test.PSFC[mytimeidx,:,:]/100
tkel=wrf.getvar(file,'tk',timeidx=mytimeidx)

cape_2d_way2=wrf.cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa,ter_follow=True)

<end code>

 

[lsun]

Thank you for the insight, and congratulations that you figured this out by yourself.

I can have a test on my end when Derecho/Casper is back online (they are under maintenance until Friday). I don't recall if I compared "wrf.getvar" with "wrf.cape_2d" in Python.

 

[epotter1]

I've been comparing different functions to calculate these values, and I suspect the reason why the most unstable cape is smaller than the surface cape might be due to the most unstable cape being calculated from the average of a 500 m parcel (at least in the wrf-python function, I think it's the same as NCL), rather than one point, and often being at or near the surface.

I'll add the reasoning here just in case anyone else is trying to compare these.

I'm comparing the wrf-python (should be the same as NCL) functions to the metpy python module, which has given some clues to this.

For the most unstable cape, the comparison is between wrf-python function wrf.getvar(file,'cape_2d') and the metpy function metpy.calc.most_unstable_cape_cin(pressure, temperature, dewpoint, height=height, depth=3000*units.m)

For the surface values, the comparison is between the lowest vertical level of wrf.getvar('file','cape_3d') and the metpy function metpy.calc.surface_based_cape_cin(pressure, temperature, dewpoint).

The wrf-python and metpy look quite similar in both cases. The surface cape values are very similar between metpy and wrf-python. The most unstable cape values follow the same patterns with a slight difference between the two. I'm putting this down to wrf-python calculating the most unstable cape from a 500 m parcel average, rather than just one point (which is the case at the surface). The cin values follow the same patterns but have a slightly bigger difference for surface cin which I can't explain.

Comparing the values of mcape with surface values, in both metpy and wrf-python the values are very close to each other, suggesting that in this case the most unstable cape is often near or at the surface. In the metpy calculation, almost all the values of maximum cape are (slightly) greater than the metpy surface cape values. I don't know why it's not all of them, but perhaps because the metpy function doesn't include the surface elevation? I'm testing over a mountainous region. In the wrf-python functions, as you say, the values of surface cape are mostly slightly higher than the values of most unstable cape. In situations where the most unstable cape is at the surface, this would be comparing a point value of most unstable cape to the average of a 500 m parcel, where the cape for other points in the parcel is lower.

This is mostly guesswork but I hope it makes some sense!