Moist Brunt-Vaisala frequency

[mwasserstein]

There is an excellent function to calculate the Brunt-Vaisala frequency in Metpy:

https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.brunt_vaisala_frequency.html

I've been searching for code to compute the moist Brunt-Vaisala frequency (a challenging calculation), and I haven't been able to find anything. Will this calculation be implemented in future versions of Metpy. Any recommendations for making this calculation myself?

[kgoebber]

We would be interested in serving a calculation such as the moist Brunt-Vaisala frequency and would welcome working with you to develop the appropriate code to go along with that. First it would help to know which version of the equation you wish to implement. A quick search yielded a paper by Durran and Klemp (1982; https://doi.org/10.1175/1520-0469(1982)039%3C2152:OTEOMO%3E2.0.CO;2) and they appear to support their equation 5.

$$N_m^{2} = \frac{g}{T}\big (\frac{dT}{dz} + \Gamma_m \big ) \big (1 + \frac{L q_s}{RT} \big ) - \frac{g}{1+q_w}\frac{dq_w}{dz}$$

Once we know which equation we want to implement, then we can open an issue and help you begin the process of using various MetPy functions and other pieces to bring it all together. In a quick inspection of the equation, it appears that this will have to be a function that works for only 1D data due to a dependence on the moist adiabatic lapse rate. So as long as you are not want to calculate this on a gridded field, you should be able to pull all of the pieces together.

[mwasserstein]

Another option that is likely easier to calculate and can probably be done on a gridded field is equation 36 from Durran and Klemp, 1982:

$$N_m^{2} = g[\frac{1 + (L q_s/ R T)}{1 + (\epsilon L^2 q_s / c_p R T^2)} \times (\frac{d \text{ln}\theta}{dz} + \frac{L}{c_p T}\frac{d q_s}{dz}) - \frac{dq_w}{dz} ]$$

The authors state that this is a "very good approximation" for the saturated Brunt-Vaisala frequency

[mwasserstein]

I've toyed around with some code to implement this on my own, but I'm not getting the results that I'd expect. Here's what I've done so far (using several metpy functions):

`

def compute_N_m_squared(ds):
# Computes N_m_squared from the ERA5 reanalysis dataset

# Read the era5 shape
num_time = len(ds.time)
num_lati = len(ds.latitude)
num_long = len(ds.longitude)
num_levl = len(ds.level)

Lv = const.Lv
R = const.Rd
eps = const.epsilon
cp = const.Cp_d
g = const.g

P = units.Quantity(np.repeat(np.repeat(np.tile(ds['level'].values[::], num_time).reshape(num_time, num_levl), num_lati), num_long).reshape(num_time, num_levl, num_lati, num_long), 'millibar')
#P = P.flatten()
T = ds.T
RH = ds.R

# Extract geopotential and assign units
geopotential =  units.Quantity(ds['Z'].values, 'meter ** 2 / second ** 2')

# Compute actual height from geopotential
height = mpcalc.geopotential_to_height(geopotential)

qs = mpcalc.saturation_mixing_ratio(total_press=P, temperature=T)
qw = mpcalc.mixing_ratio_from_relative_humidity(pressure = P, temperature = T, relative_humidity=RH)
theta = mpcalc.potential_temperature(P, T)

# Term one
numerator = 1 + ( (Lv * qs) / (R * T * units.kelvin) )
denominator = 1 + ( (eps * (Lv ** 2) * qs) / (cp * R * (T ** 2) * (units.kelvin ** 2)) )

term_1 = numerator / denominator

# Term two
# Term 2_a
term_2_a = mpcalc.first_derivative(np.log(theta.values), x = height)

# Term 2_b_1
term_2_b1 = Lv / (cp * T * units.kelvin)
# Term 2_b_2
term_2_b2 = mpcalc.first_derivative(qs, x = height)

term_2 = term_2_a + (term_2_b1 * term_2_b2)

# Term 3
term_3 = mpcalc.first_derivative(qw, x = height)

N_m_squared = g * (term_1 * term_2 - term_3 )

return N_m_squared

`

[kgoebber]

Okay, so I think there are some things going on with multiplying values (e.g., constants) against a DataArray that was cause some unit issues. I have a full working example with GFS data...

from datetime import datetime

import metpy.constants as mpconst
import xarray as xr

def compute_N_m_squared(P, T, RH, GEOPOT):
    import metpy.calc as mpcalc
    import metpy.constants as mpconst
    from metpy.units import units
    # Computes N_m_squared from the ERA5 reanalysis dataset

    Lv = mpconst.Lv
    R = mpconst.Rd
    eps = mpconst.epsilon
    cp = mpconst.Cp_d
    g = mpconst.g
    
    # Compute actual height from geopotential
    height = mpcalc.geopotential_to_height(GEOPOT)
    
    qs = mpcalc.saturation_mixing_ratio(total_press=P.data[:, None, None] * units(P.units), temperature=T)
    qw = mpcalc.mixing_ratio_from_relative_humidity(pressure = P.data[:, None, None] * units(P.units),
                                                    temperature = T, relative_humidity=RH)
    theta = mpcalc.potential_temperature(P.data[:, None, None] * units(P.units), T)
    
    # Term one
    numerator = 1 + ( (Lv * qs) / (R * T) )
    denominator = 1 + ( (eps * (Lv ** 2) * qs) / (cp * R * (T ** 2)) )

    term_1 = numerator / denominator

    # Term two
    # Term 2_a
    term_2_a = mpcalc.first_derivative(theta, x = height) / theta

    # Term 2_b_1
    term_2_b1 = Lv / (cp * T)
    # Term 2_b_2
    term_2_b2 = mpcalc.first_derivative(qs, x = height)

    term_2 = term_2_a + (term_2_b1 * term_2_b2)

    # Term 3
    term_3 = mpcalc.first_derivative(qw, x = height)

    N_m_squared = g * (term_1 * term_2 - term_3 )

    return N_m_squared

date = datetime(2023, 1, 15, 12)

# Get GFS data for contouring
ds = xr.open_dataset('https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/'
                     f'Global_onedeg_ana/GFS_Global_onedeg_ana_{date:%Y%m%d_%H%M}.grib2')

# Subset data to be just over the CONUS
ds = ds.metpy.sel(time=date)

# Temperature, RH, and Geopotential Height as Unit Arrays
temperature = ds.Temperature_isobaric.metpy.quantify()
relative_humidity = ds.Relative_humidity_isobaric.metpy.quantify()
geopotential = ds.Geopotential_height_isobaric.metpy.quantify() * mpconst.g

# Pressure as a DataArray
pressure = ds.Temperature_isobaric.metpy.vertical

compute_N_m_squared(pressure, temperature, relative_humidity, geopotential)*1e4

Try to modify with ERA5 variable names and/or see if this gives sensible values from the GFS output.

[kgoebber]

Okay, I wonder if we've got something going on with reversal of the pressure coordinate. The following code uses GFS data from the RDA archive at NCAR appears to work well. Only one modification (I think) from prior code and that is to have qw = qs + ql which seems to align better with the description from the reference paper.

from datetime import datetime

import metpy.constants as mpconst
from metpy.plots import SkewT
import xarray as xr

def compute_N_m_squared(P, T, RH, GEOPOT):
    import numpy as np
    import metpy.calc as mpcalc
    import metpy.constants as mpconst
    from metpy.units import units
    # Computes N_m_squared from the ERA5 reanalysis dataset

    Lv = mpconst.Lv
    R = mpconst.Rd
    eps = mpconst.epsilon
    cp = mpconst.Cp_d
    g = mpconst.g
    
    # Compute actual height from geopotential
    height = mpcalc.geopotential_to_height(GEOPOT)

    qs = mpcalc.saturation_mixing_ratio(total_press=P.values[:, None, None] * units(P.units), temperature=T)
    ql = mpcalc.mixing_ratio_from_relative_humidity(pressure = P.values[:, None, None] * units(P.units),
                                                    temperature = T, relative_humidity=RH)
    qw = qs + ql

    theta = mpcalc.potential_temperature(P.values[:, None, None] * units(P.units), T)
    
    # Term one
    numerator = 1 + ( (Lv * qs) / (R * T) )
    denominator = 1 + ( (eps * (Lv ** 2) * qs) / (cp * R * (T ** 2)) )

    term_1 = numerator / denominator

    # Term two
    # Term 2_a
    term_2_a = mpcalc.first_derivative(theta, x = height) / theta

    # Term 2_b_1
    term_2_b1 = Lv / (cp * T)
    # Term 2_b_2
    term_2_b2 = mpcalc.first_derivative(qs, x = height)

    term_2 = term_2_a + (term_2_b1 * term_2_b2)

    # Term 3
    term_3 = mpcalc.first_derivative(qw, x = height)

    N_m_squared = g * (term_1 * term_2 - term_3 )

    return N_m_squared

date = datetime(2005, 1, 10, 18)

# Get GFS data
ds = xr.open_dataset('https://rda.ucar.edu/thredds/dodsC/aggregations/g/ds083.2/1/TP')

subset = dict(time=date, vertical=np.arange(100, 1000, 50)*units.hPa)

# Temperature, RH, and Geopotential Height as Unit Arrays
temperature = ds.Temperature_isobaric.metpy.sel(subset).metpy.quantify()

relative_humidity = ds.Relative_humidity_isobaric.metpy.sel(subset).metpy.quantify()
# Rename a vertical coordinate to facilitate calculation
relative_humidity = relative_humidity.rename({relative_humidity.metpy.vertical.name:temperature.metpy.vertical.name})
dewpoint = mpcalc.dewpoint_from_relative_humidity(temperature, relative_humidity)

geopotential = ds.Geopotential_height_isobaric.metpy.sel(subset).metpy.quantify() * mpconst.g

uwnd = ds['u-component_of_wind_isobaric'].metpy.sel(subset)
vwnd = ds['v-component_of_wind_isobaric'].metpy.sel(subset)

# Pressure as a DataArray
pressure = ds.Temperature_isobaric.metpy.sel(subset).metpy.vertical

# Compute Moist Brunt-Vaisala Frequency
nm2 = compute_N_m_squared(pressure, temperature, relative_humidity, geopotential)*1e4

# Compute the Brunt Vaisala Frequency from MetPy
n2 = mpcalc.brunt_vaisala_frequency_squared(mpcalc.geopotential_to_height(geopotential),
                                            mpcalc.potential_temperature(pressure, temperature))*1e4

# Make Plot
point_ilat = 50
point_ilon = 50

fig = plt.figure(figsize=(15, 10))
ax = plt.subplot(121)
ax.invert_yaxis()
ax.semilogy(nm2[:, point_ilat, point_ilon], pressure.metpy.convert_units('hPa'), label=r'$N_m^2$')
ax.semilogy(n2[:, point_ilat, point_ilon], pressure.metpy.convert_units('hPa'), label=r'$N^2$')
plt.legend()
ax.set_ylabel('Pressure (hPa)')
ax.set_yticks(np.arange(1000, 99, -100), np.arange(1000, 99, -100))
ax.set_ylim(1050, 100)
ax.set_xlabel(r'N^2*1e4')

ax.set_title(f'BV vs. Moist BV Calculation {date}')

skew = SkewT(fig, rotation=45, subplot=(1,2,2))
skew.plot(pressure.metpy.convert_units('hPa'),
          temperature[:, point_ilat, point_ilon].metpy.convert_units('degC'), color='red')
skew.plot(pressure.metpy.convert_units('hPa'),
          dewpoint[:, point_ilat, point_ilon].metpy.convert_units('degC'), color='green')
skew.plot_barbs(pressure.metpy.convert_units('hPa'),
                uwnd[:, point_ilat, point_ilon].metpy.convert_units('knot'),
                vwnd[:, point_ilat, point_ilon].metpy.convert_units('knot'))
skew.plot_moist_adiabats()
plt.xlim(-30, 20)
plt.show()

What particular location are from your examples?

[mwasserstein]

Hi Kevin. Thanks again for all your help. I am plotting using the ERA5 data for just east of Salt Lake City, UT. (40.5, -111.75). I have tried the computation using your modifications, and it still looks off. Here is an example. The first shows the GFS for the gridpoint (40, -112) and the second shows the ERA5 for the gridpoint (40.5, -111.75):

I think that part of what could be causing this issue is that the ERA5 has RH values > 100% (and > 120% in some cases!). We are using equation 36 from Durran and Klemp, 1982, which is "very good approximation" for the saturated Brunt-Vaisala frequency. Maybe this approximation just does not do well when RH > 100%?

[kgoebber]

Interesting. RH > 100% could be it. You could force the ERA5 to have RH <= 100 with the following code:

relative_humidity.values[(relative_humidity > 100 * units.percent).values] = 100 * units.percent

To see if that changes the behavior at all.

[DanielAdriaansen]

Hello, I don't have much to add regarding the science here, but the location of your grid point reminded me of this:
https://confluence.ecmwf.int/pages/viewpage.action?pageId=153398938'

It looks like the lowest pressure level in the GFS data is higher up than the ERA5 data, which appear to go all the way to 1000 hPa. Perhaps the ERA5 below surface data is causing an issue (which may lead to RH >> 100%)?

[mwasserstein]

Kevin, I tried forcing the RH to <= 100, and that did not solve my issue. I'm really not sure what is causing these strange values. And Daniel, that is a good thought. In the grid point I am using, there is subterranean data, which could be screwing things up. Maybe the ERA5 is just not the best dataset for my application.

I do think the code that we have to calculate $N_m^2$ is correct, though.