I'm very new to metpy. As a first example I tried to calculate the LCL of an air parcel with T=25°C, p=950hPa and relative humidity=70%.
First, the mixing ratio of this air is calculated to be
mr: 15.016493919350289 gram / kilogram
Afterwards, the LCL is calculated as:
pLCL: 871.9597894290378 hectopascal
TLCL: 17.786684840244334 degree_Celsius
Now, I tried to calculate saturation equivalent potential temperature and mixing_ratio for levels starting with pLCL. To my surprise, the mixing ratio, calculated with mixing_ratio_from_relative_humidity
at the LCL differs from the mixing ratio of the original parcel. Instead of getting again mr=15.016493919350289 gram / kilogram I get 14.862703501798386 gram / kilogram. Shouldn't it be the same, because at LCL the relative humidity is 100% and during dry adiabatic rise to LCL mixing ratio cannot change?
There is a link to this question in: https://earthscience.stackexchange.com/questions/25188/why-is-mixing-ratio-at-lcl-not-the-same-as-for-starting-condition
from metpy.calc import dewpoint_from_relative_humidity
from metpy.calc import equivalent_potential_temperature
from metpy.calc import mixing_ratio_from_relative_humidity
from metpy.calc import equivalent_potential_temperature
from metpy.calc import lcl
from metpy.calc import moist_lapse
from metpy.calc import saturation_equivalent_potential_temperature
from metpy.units import units
p = 950 * units.hPa
T = 25 * units.degC
rh = 70 * units.percent
print('p:\t', p)
print('T:\t', T)
print('RF:\t', rh)
dp = dewpoint_from_relative_humidity(T, rh)
mr = mixing_ratio_from_relative_humidity(p, T, rh).to('g/kg')
print('mr:\t', mr)
ept = equivalent_potential_temperature(p, T, dp)
print('Θe :\t', ept.to(units.degC))
mylcl = lcl(p, T, dp)
plcl = mylcl[0]
tlcl = mylcl[1]
print('pLCL:\t', plcl)
print('TLCL:\t', tlcl)
plevs = [plcl.magnitude, 800, 700, 600, 500, 400, 300, 200, 100, 50, 25] * units.hPa
ml = moist_lapse(plevs, tlcl).to('degC')
print()
print('Moist adiabatic lift, starting from LCL')
print()
rh100 = 100 * units.percent
for idx in range(len(plevs)):
p = plevs[idx]
T = ml[idx]
print('-------------------------------------------------------------')
print('p:\t', p)
print('T(p):\t', T)
ept = saturation_equivalent_potential_temperature(plevs[idx], ml[idx])
print('Θe(p):\t', ept.to(units.degC))
mr = mixing_ratio_from_relative_humidity(p, T, rh100)
print('mr:\t', mr.to('g/kg'))
p: 950 hectopascal
T: 25 degree_Celsius
RF: 70 percent
mr: 15.016493919350289 gram / kilogram
Θe : 73.59933073452271 degree_Celsius
pLCL: 871.9597894290378 hectopascal
TLCL: 17.786684840244334 degree_Celsius
Moist adiabatic lift, starting from LCL
-------------------------------------------------------------
p: 871.9597894290378 hectopascal
T(p): 17.786684840244334 degree_Celsius
Θe(p): 73.56266438016314 degree_Celsius
mr: 14.862703501798386 gram / kilogram <--- !!!!!!!
-------------------------------------------------------------
p: 800.0 hectopascal
T(p): 14.680478970030833 degree_Celsius
Θe(p): 73.7567612561337 degree_Celsius
mr: 13.254531670351966 gram / kilogram
-------------------------------------------------------------
p: 700.0 hectopascal
T(p): 9.70674918909208 degree_Celsius
Θe(p): 73.97100713514288 degree_Celsius
mr: 10.878278101841651 gram / kilogram
-------------------------------------------------------------
p: 600.0 hectopascal
T(p): 3.6628633969115754 degree_Celsius
Θe(p): 74.0966050124087 degree_Celsius
mr: 8.34266379949833 gram / kilogram
-------------------------------------------------------------
p: 500.0 hectopascal
T(p): -4.036747412608861 degree_Celsius
Θe(p): 74.09741860835788 degree_Celsius
mr: 5.6959850509157075 gram / kilogram
-------------------------------------------------------------
p: 400.0 hectopascal
T(p): -14.539107430252784 degree_Celsius
Θe(p): 73.9329618099161 degree_Celsius
mr: 3.10991709553994 gram / kilogram
-------------------------------------------------------------
p: 300.0 hectopascal
T(p): -30.169336134705304 degree_Celsius
Θe(p): 73.61031144188621 degree_Celsius
mr: 1.043016676418755 gram / kilogram
-------------------------------------------------------------
p: 200.0 hectopascal
T(p): -54.703600518566134 degree_Celsius
Θe(p): 73.32450117793161 degree_Celsius
mr: 0.11362204759565221 gram / kilogram
-------------------------------------------------------------
p: 100.0 hectopascal
T(p): -93.72997677325506 degree_Celsius
Θe(p): 73.25906799180353 degree_Celsius
mr: 0.0005989154996261802 gram / kilogram
-------------------------------------------------------------
p: 50.0 hectopascal
T(p): -125.96434978295119 degree_Celsius
Θe(p): 73.2583328918671 degree_Celsius
mr: 4.5360156886986637e-07 gram / kilogram
-------------------------------------------------------------
p: 25.0 hectopascal
T(p): -152.4084163524267 degree_Celsius
Θe(p): 73.25830368606137 degree_Celsius
mr: 2.1998938295740055e-11 gram / kilogram
You didn't share how you calculated the second mixing ratio value, but I assume it's mixing_ratio_from_relative_humidity(plcl, tlcl, 100 * units.percent)
. Regardless, I think this stems from an inconsistency in MetPy in how relative humidity (and thus saturation) is assumed, as noted in this GitHub issue. The current implementation of mixing_ratio_from_relative_humidity()
relies upon the WMO definition (since it's more direct for this calculation), which is: RH = w / w_s
(where w
is mixing ratio and _s
denotes the quantity at saturation.
lcl()
on the other hand, assumes saturation is reached when the ambient water vapor pressure reaches its saturation value. So in this case, effectively using a definition of RH = e / e_s
(where e
is water vapor partial pressure)--the definition given by the AMS glossary and used in a lot of other software, like GEMPAK.
You can see the conservation of mixing ratio if you change your initial value to go through a computation path that uses the AMS definition of RH:
Td = mpcalc.dewpoint_from_relative_humidity(T, rh)
mr = mpcalc.mixing_ratio(mpcalc.saturation_vapor_pressure(Td), p)
which gives the 14.86 g/kg value, and also then matches the mixing ratio at the LCL:
mpcalc.mixing_ratio(mpcalc.saturation_vapor_pressure(tlcl), plcl)
At some point we'll resolve the GitHub issue and make things internally consistent to avoid confusion like this. It should be noted, though, that your original values are within ~1% of each other, so we're not talking about a significant discrepancy.