Lab 2-3: Atmospheric Thermodynamics and Stability

SAIL Radiosonde Data

This notebook is designed to investigate 2022 winter radiosonde data from SAIL launched from Gothic, CO over a variety of conditions

The act package was applied for this notebook to download and easily visualize the data. While not required for class, you can check out the documentation here.

# Import packages to be used
import os

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

%matplotlib inline

Examine a sounding

Let’s start by taking a look at a sounding from a clear, calm day on 10 January 2022. You can check how billy barr described the day’s weather by looking here. Note that Gothic, Colorado is an old mining town in the mountains of Colorado at 2,900 m (9,500 ft) elevation. Note that the code below expects to find the sounding data inside folder labeled “data” in your working folder. If you haven’t already, make a folder named data and put the .csv files in that folder.

# first, read in the data from the sounding

# Get path to data
# filepath= os.path.join(os.getcwd(),'data') #use this if you put the data in a data folder
# use the below if you put the file in the same folder as this notebook:
filepath= os.path.join(os.getcwd())
date = '2022-01-10'

# Read the data into a file
sonde_temp_profile = pd.read_csv(os.path.join(filepath,f'{date}_radiosonde.csv'))

# take a look at what is included in the file
sonde_temp_profile.head()

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[3], line 10
      6 filepath= os.path.join(os.getcwd())
      7 date = '2022-01-10'
      8 
      9 # Read the data into a file
---> 10 sonde_temp_profile = pd.read_csv(os.path.join(filepath,f'{date}_radiosonde.csv'))
     11 
     12 # take a look at what is included in the file
     13 sonde_temp_profile.head()

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/pandas/io/parsers/readers.py:873, in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, skip_blank_lines, parse_dates, date_format, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, on_bad_lines, low_memory, memory_map, float_precision, storage_options, dtype_backend)
    861 kwds_defaults = _refine_defaults_read(
    862     dialect,
    863     delimiter,
   (...)    869     dtype_backend=dtype_backend,
    870 )
    871 kwds.update(kwds_defaults)
--> 873 return _read(filepath_or_buffer, kwds)

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/pandas/io/parsers/readers.py:300, in _read(filepath_or_buffer, kwds)
    297 _validate_names(kwds.get("names", None))
    299 # Create the parser.
--> 300 parser = TextFileReader(filepath_or_buffer, **kwds)
    302 if chunksize or iterator:
    303     return parser

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/pandas/io/parsers/readers.py:1645, in TextFileReader.__init__(self, f, engine, **kwds)
   1642     self.options["has_index_names"] = kwds["has_index_names"]
   1644 self.handles: IOHandles | None = None
-> 1645 self._engine = self._make_engine(f, self.engine)

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/pandas/io/parsers/readers.py:1904, in TextFileReader._make_engine(self, f, engine)
   1902     if "b" not in mode:
   1903         mode += "b"
-> 1904 self.handles = get_handle(
   1905     f,
   1906     mode,
   1907     encoding=self.options.get("encoding", None),
   1908     compression=self.options.get("compression", None),
   1909     memory_map=self.options.get("memory_map", False),
   1910     is_text=is_text,
   1911     errors=self.options.get("encoding_errors", "strict"),
   1912     storage_options=self.options.get("storage_options", None),
   1913 )
   1914 assert self.handles is not None
   1915 f = self.handles.handle

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/pandas/io/common.py:926, in get_handle(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)
    921 elif isinstance(handle, str):
    922     # Check whether the filename is to be opened in binary mode.
    923     # Binary mode does not support 'encoding' and 'newline'.
    924     if ioargs.encoding and "b" not in ioargs.mode:
    925         # Encoding
--> 926         handle = open(
    927             handle,
    928             ioargs.mode,
    929             encoding=ioargs.encoding,
    930             errors=errors,
    931             newline="",
    932         )
    933     else:
    934         # Binary mode
    935         handle = open(handle, ioargs.mode)

FileNotFoundError: [Errno 2] No such file or directory: '/home/runner/work/Fluid_Flows/Fluid_Flows/book/modules/module2/2022-01-10_radiosonde.csv'

We can see that the sonde was launched at 11:34, and its first recording was taken at an atmospheric pressure of 714 hPa. We can also look at specific values. Let’s pretend that the first value represents the surface air (it’s actually slightly above the surface, but it’s close). Then can we calculate the characteristics of the surface parcel of air.

Tsurf = sonde_temp_profile.loc[0,['tdry']]
RHsurf = sonde_temp_profile.loc[0,['rh']]
print(Tsurf)
print(RHsurf)
# You should see temperature in deg C and RH in percent
# Also note that we get information about the data type (dtype) because this is a pandas dataframe

tdry   -9.97983
Name: 0, dtype: object
rh    52.05915
Name: 0, dtype: object

Saturated vapor pressure is a function of air temperature. Recall that the formula, as shown here, calculates saturated vapor pressure in units of kPa with an input of temperature in degrees C. $\( e_{sat} = 0.6108*exp(\frac{17.27T}{237.3+T}) \)$

# We can use our function for the saturated vapor pressure and the python math forumulas 
# to calculate the saturated vapor pressure esat near our surface
# The first step is to change the datatype of the Tsurf variable from above so that it is a floating point
# number and can be used in equations.

T = np.asarray(Tsurf).astype(np.float64)
print(T)

[-9.97983]

T1 = 10
T2 = 2
esat1 = 0.6108*np.exp((17.27*T)/(237.3+T))
# this formula outputs esat in kPa
print(esat1)

[0.28616842]

We can also use this formula to create and plot the saturated vapor pressure curve.

T_curve=list(range(-40,40))
T_curve=np.asarray(T_curve)
# We turn it into a numpy array so that we can use numpy to calculate esat for the entire array
print(T_curve)

[-40 -39 -38 -37 -36 -35 -34 -33 -32 -31 -30 -29 -28 -27 -26 -25 -24 -23
 -22 -21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10  -9  -8  -7  -6  -5
  -4  -3  -2  -1   0   1   2   3   4   5   6   7   8   9  10  11  12  13
  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31
  32  33  34  35  36  37  38  39]

# Now, we use numpy with the equation over our entire array
esat_curve = 0.6108*np.exp((17.27*T_curve)/(237.3+T_curve))

# Create a new figure.
plt.figure()

# Use the plot() function to plot the year on the x-axis, peak flow values on
# the y-axis with an open circle representing each peak flow value.
plt.plot(T_curve, # our x value
         esat_curve # our y value
         )
# Label the axes and title.
plt.xlabel('Temperature (deg C)')
plt.ylabel('Saturated Vapor Pressure (kPa)')
plt.title('(Fig. 1) Saturated Vapor Pressure Curve');

We know that Relative Humidity, $\( RH = e_v/e_{sat}*100% \)\( where \)e_v\( is vapor pressure and \)e_{sat}$ is saturated vapor pressure.

# We can combine the measured RH and T to determine the actual vapor pressure (in kPa)
RH = np.asarray(RHsurf).astype(np.float64)
ev = esat1*RH/100
# Note that this will be in the same units as esat, which in this case is kPa
print(ev)

[0.14897685]

With this information, as well as the measurement of total air pressure, we can calculate all of the other moisture-related variables.
Dew point temperature, in deg C can be defined as, $\( T_d = \frac{ln(e_v) + 0.49299}{0.0707-0.00421ln(e_v)} \)$

and can also be defined by the Magnus formulation (Murray 1967) as, in deg C,

\[ T_d = \frac{243.04[ln(\frac{RH}{100}) + \frac{17.625T}{243.04+T}}{17.625 - ln(\frac{RH}{100})-\frac{17.625T}{243.04+T}} \]

#We can use our measurements from above in these calculations and see how well they match.
Td1 = (np.log(ev)+0.49299)/(0.0707-0.00421*np.log(ev))
print(Td1)

[-17.92494446]

Td2 = (243.04*(np.log(RH/100)+(17.625*T)/(243.04+T)))/(17.625 - np.log(RH/100) - (17.625*T)/(243.04+T))
print(Td2)

[-17.9734832]

If we also look at the total pressure measured by the radiosonde, we can calculate the mixing ratio $\( R = \frac{0.622e_v}{p-e_v} \)$ Remember that the mixing ratio is the ratio of the mass of water vapor to the mass of dry air only.

# First, convert the bottom pressure level from our data frame into a floating point number
Psurf = sonde_temp_profile.loc[0,['pres']]
P = np.asarray(Psurf).astype(np.float64)
# Check that we have these in the same units.  P from the sonde is in mb or hPa and ev is in kPa
# if we look online, we see that 1 kPa = 10 mb, to keep everything in kPa, we need to divide P by 10
# now everything is in kPa
R = 0.622*ev/((P/10)-ev)
print(R)

[0.00130091]

In class, we talked about the mixing ratio being approximately equal to the specific humidity and being approximately equal to 0.622(ev/p) . We can check if this is true.

R2 = 0.622*ev/((P/10))
print(R2)

[0.00129819]

The specific humidity, q is defined as the ratio of the mass of vapor to the mass of the total air (dry plus wet), as $\( q = \frac{\rho_v}{\rho_a} = \frac{\rho_v}{(\rho_v+\rho_d)} \)\( where \)\rho_a\( is the density of the total air, \)\rho_v\( is the denisty of the water vapor alone and \)\rho_d$ is the density of dry air alone.

The radiosonde doesn’t directly measure density, so to fully calculate this, we need to apply the ideal gas law and Dalton’s law of partial pressures. Following from page 17 of Shuttleworth, we get $\( q = \frac{0.62e_v}{p-0.38e_v} \)$

q = 0.62*ev/((P/10)-0.38*ev)
print(q)

[0.00129504]

As stated in the textbook and in class, these are all very similar to each other.

Before we consider buoyancy, stability, and vertical motion, we need one more term. This is necessary because while warmer air is less dense than cooler air, moister air is less dense than drier air. Virtual temperature: $\( T_v = T(1+0.61q) \)\( is the temperature that dry air would have if it had the same density and temperature as the moist air, where in this formula, \)T\( and \)T_v$ are in Kelvin. It allows us to consider moisture effects on density.

#Let's calculate the virtual temperature of our surface air from this sounding.
# Note that this equation applies to Kelvin, which is 273.15 greater than deg C
TvK = (T+273.15)*(1+0.61*q)
print(TvK)
print(TvK-273.14)

[263.37806823]
[-9.76193177]

Profiles and Vertical Motion

Let’s continue with this same sounding to explore vertical profiles and tendencies for vertical motion.

Buoyant acceleration is defined as $\( a_p = -g\frac{(\rho_a' - \rho_a)}{\rho_a} \)$

or

\[ a_p = -g\frac{(T_v - T_v')}{T_v} \]

where the ‘ indicates the perturbed parcel, and the non-prime value indicates the surrounding air. Here we see that a lower density will move up, and a warmer temperature will move up.

# for our January 10th sounding, let's plot the actual temperature, 
# the virtual temperature, and the temperature of a parcel following the dry adiabatic lapse rate.

# Create a new figure.
# Start with the actual temperature, which we measured with the sounding
fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(7,7), tight_layout=True)
axs.plot(sonde_temp_profile['tdry'],sonde_temp_profile['alt'])
axs.set_ylim(3000,6000)
axs.set_xlim(-20,0)
axs.set_xlabel('temperature (C)')
axs.set_ylabel('altitude (m)')
plt.title('(Fig. 2) Atmospheric Temperature Profile from 10 January 2022');

# calculate the dry adiabat with elevation
# first determine the surface elevation
Altsurf = sonde_temp_profile.loc[0,['alt']]
Altsurf = np.asarray(Altsurf).astype(np.float64)
print(Altsurf)
print(T)
alt4adiabat =np.asarray(list(range(3055,6055))) 
T4adiabat = T - 9.8*(alt4adiabat-3055)/1000

# And add this to the figure above
fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(7,7), tight_layout=True)
axs.plot(sonde_temp_profile['tdry'],sonde_temp_profile['alt'])
axs.plot(T4adiabat,alt4adiabat,'r-')
axs.set_ylim(3000,6000)
axs.set_xlim(-20,0)
axs.set_xlabel('temperature (C)')
axs.set_ylabel('altitude (m)')
plt.title('(Fig. 2) Atmospheric Temperature Profile & Dry Adiabat from 10 Jan 2022');

[3055.1375]
[-9.97983]

Do we have stable or unstable conditions in this sounding?

# Now we use the data in the sounding to calculate the virtual temperature at each altitude
# Note that what we are doing with the code is using the equations presented above but iterating over all
# of the values at all of the heights in the sounding dataframe, so we get arrays of numbers
esat = 0.6108*np.exp((17.27*sonde_temp_profile['tdry'])/(237.3+sonde_temp_profile['tdry']))
ev = esat*sonde_temp_profile['rh']/100
q = 0.62*ev/((sonde_temp_profile['pres']/10)-0.38*ev)
TvK = (sonde_temp_profile['tdry']+273.15)*(1+0.61*q)
#TvK = (T+273.15)*(1+0.61*q)
print(TvK)

#Now, think about how you would modify the code presented to add it to the plot yourself.
# You can use the + button to put a new cell below and then copy and past the code from above into that cell.
# Think about the units of the plot above.  Will TvK plot on the same graph as is?
# Think about the vertical coordinate for TvK -- is it the same as for the adiabat or for tdry?

0     263.378068
1     267.845751
2     268.835092
3     268.248066
4     267.485153
         ...    
91    216.796763
92    217.033717
93    218.668511
94    220.611648
95    221.227319
Length: 96, dtype: float64

Disclaimer, solution is presented below. Please try to modify the code yourself before you look.

# 

#Now, modify the code presented to add it to the plot yourself.
# And add this to the figure above
fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(7,7), tight_layout=True)
axs.plot(sonde_temp_profile['tdry'],sonde_temp_profile['alt'])
axs.plot(T4adiabat,alt4adiabat,'r-')
axs.plot(TvK-273.15,sonde_temp_profile['alt'],'k-')
axs.set_ylim(3000,6000)
axs.set_xlim(-20,0)
axs.set_xlabel('temperature (C)')
axs.set_ylabel('altitude (m)')
plt.title('(Fig. 3) T Profile (blue), Tvprofile (black), & Dry Adiabat (red) from 10 Jan 2022');

So, what difference did Tvirtual make?

First, the difference is slight. Is there any change in the stability of this sounding? Second, considering the moisture content makes the equivalent dry air temperature a bit warmer. Does this make sense with your understanding of the virtual temperature?
Given the moist air is less dense than dry air, and warm air is less dense than cold air, more moisture would make the virtual temperature warmer than the measured temperature.

Plotting function and examples

We have included four days of sounding data. To make plotting easier, we can define a function that takes inputs and creates a plot. You may want to use this for your homework.

def plot_sonde(filepath,date,variable, limits, title, unit='$\degree$C'):
    # Apply date filter
    sonde_temp_profile = pd.read_csv(os.path.join(filepath,f'{date}_radiosonde.csv'))
        
    # Plot data
    fig, axs = plt.subplots(ncols=2, figsize=(10,5), tight_layout=True)
    axs[0].plot(sonde_temp_profile[variable],sonde_temp_profile['alt'])
    axs[0].set_xlim(limits[0],limits[1])
    axs[0].set_ylim(3000,10000)
    axs[0].set_xlabel(f'{variable.upper()} ({unit})');
    axs[0].set_ylabel('Altitude (m)');
    axs[0].set_title(f'{variable.upper()} with Altitude');

    axs[1].plot(sonde_temp_profile[variable],sonde_temp_profile['pres'])
    axs[1].set_xlim(limits[0],limits[1])
    axs[1].set_ylim(725,350)
    axs[1].set_xlabel(f'{variable.upper()} ({unit})');
    axs[1].set_ylabel('Pressure (hPa)');
    axs[1].set_title(f'{variable.upper()} with Pressure');

    fig.suptitle(f'{title} on {date}')
    return

<>:1: SyntaxWarning: invalid escape sequence '\d'
<>:1: SyntaxWarning: invalid escape sequence '\d'
/tmp/ipykernel_128/2951637969.py:1: SyntaxWarning: invalid escape sequence '\d'
  def plot_sonde(filepath,date,variable, limits, title, unit='$\degree$C'):

We can use this function to generate plots comparing temperature with height on different days, under different conditions:

# Get path to data:  USE THIS if you have your radiosonde files in a data folder
#filepath= os.path.join(os.getcwd(),'data')
# Use the code below if you just copied the data files into the same folder as this notebook
filepath= os.path.join(os.getcwd())

# Apply date filter
date = '2022-01-10'
plot_sonde(filepath, date, 'tdry', (-50,0), title='Clear Conditions with Inversion')

# Apply date filter
date = '2022-01-03'
plot_sonde(filepath, date, 'tdry', (-50,0), title='Clear Conditions with Inversion')

# Apply date filter
date = '2021-12-26'
plot_sonde(filepath, date, 'tdry', (-50,0), title='Stormy/Cloudy Conditions without Inversion')

# Apply date filter
date = '2022-01-06'
plot_sonde(filepath, date, 'tdry', (-50,0), title='Stormy/Cloudy Conditions without Inversion')

Plots comparing relative humidity with height during different conditions

# Apply date filter
date = '2022-01-06'
plot_sonde(filepath, date, 'rh', (0,100), title='Stormy/Cloudy Conditions without Inversion', unit='%')

# Apply date filter
date = '2021-12-26'
plot_sonde(filepath, date, 'rh', (0,100), title='Stormy/Cloudy Conditions without Inversion', unit='%')

# Apply date filter
date = '2022-01-10'
plot_sonde(filepath, date, 'rh', (0,100), title='Clear Conditions with Inversion', unit='%')

In addition to the balloon soundings, the SAIL campaign has also placed a ceilometer at Gothic. The ceilometer uses a laser beam to measure backscatter from particles in the atmosphere. Most particles are in the lowest layer, or boundary layer, of the atmosphere, so you can detect the height of the boundary layer as the point where there are fewer returns. 100% reflection occurs at the bottom of a cloud, so a cloud layer can be identified as the height of strong returns with no returns above it. Rain or snow below a cloud layer will also cause strong returns. Below you can see ceilometer plots showing boundary layer height and cloud height with time for our study site. Look for the dates corresponding to the radiosonde data you have plotted above to check if you have identified the cloud base height correctly. (Note that ceilometer is plotting height above the surface, so you will want to add 3000 m to it to compare with the graphs above.)