Notebook 1.1: Atmospheric Modeling

January 14, 2010

Famous right-wing radio personality Rush Limbaugh claims that chlorofluorocarbons can’t possibly deplete ozone in the stratosphere, because they’re heavier than air. They’ll settle close to the ground, says Rush, where they can’t do any harm.

The Euler notebook posted below asks students to apply the barometric equation to see whether or not gases really stratify in the atmosphere based on their molar masses. To use it, download the Euler Math Toolbox. After installing the toolbox, just cut and paste the following code into a plain text file, save it with an .en file extension, and double-click the file. It should open up in the Toolbox.

I used this in my Fall ‘09 PChem class (the second in the series). Italics mark the items students had to fill in themselves. Some of my students struggled with this assignment because they couldn’t fathom doing calculations with arrays. I really should have had a notebook before this one that focused exclusively on arithmetic with arrays.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition.

---------------------------------snip here-----------------------------
Notebook 1.1 - Atmosphere Modeling
% As before, hit return on each line to execute the calculation on that
% line. Use ALT-Insert to insert a new line, and ALT-Delete to remove a
% line.
% -----------------------------------------------------------------------
% In this notebook, we'll compare two models for describing the
% change in air pressure with altitude. Both models are based on the
% barometric formula
% P = Po*exp(-h*M*g/(R*T))
% where P is the pressure, Po is the pressure at altitude 0, h is the
% altitude, M is the molar mass, g is the acceleration due to gravity,
% R is the ideal gas constant, and T is the temperature in kelvins.
% * MODEL A assumes that air behaves like an ideal gas with a
% single molar mass, Mo.
% * MODEL B assumes that each gas in air behaves as an ideal
% gas, and each gas develops its own separate altitude
% profile. In this model, we should see that air is richer in
% heavy gases at lower altitudes.
% TASK 2A. Let's begin with model A by using the barometric formula to plot
% air pressure (in atm, on the x-axis) against altitude (on the y-axis).
% First, we need to know how temperature varies with altitude.
% We can put the temperature vs. altitude data into lists called
% 'arrays'. Arrays can be entered into Euler notebooks by using
% square brackets. For example, here are the data arrays for altitude
% (h, in meters) and temperature (t, in degrees Celsius).


% Each altitude in h corresponds to a temperature in t.
% We can plot arrays and show the data as separate
% points using plot2d's 'points=1' option. We can choose the style of the
% points using the 'style' option. Style can be set to "o", ".", "x",
% or "[]" with points=1.
% Always title your plots, and label the axes. Include units in the
% axis labels.
>title("Variation of Temperature with Altitude");
>xlabel("temperature (degrees Celsius)");
>ylabel("altitude (m)");

% Next, we need the average molar mass of air at sea level, Mo (in
% kg/mol). The molecule fractions in air are 78.09% N2, 20.95%
% O2, 0.93% Ar, and 0.03% CO2 at sea level, with negligible
% amounts of other gases in dry air. We can use arrays again, with
% each separate gas as an 'element' in the array (N2 is element 1,
% O2 is element 2, Ar is element 3, and CO2 is element 4):
>percent = [78.09, 20.95, 0.93, 0.03];
>abundance = percent/100;
>mass = [2*14.0067, 2*15.9994, 39.948, 12.0107+2*15.9994];
% If we're going to do our calculations in SI units, we need the
% masses in kg/mol, not g/mol. Note that we can do the conversion
% on all elements of the array at once.
>mass = mass/1000;
% We can access the i-th element of array x by writing x[i].
% For example, the mass of N2 is mass[1]. The abundance of Ar is
% abundance[3].
>Mo = abundance[1]*mass[1] + abundance[2]*mass[2] + abundance[3]*mass[3] + abundance[4]*mass[4]
% It's much easier to use the sum function, which adds up all of the
% elements of an array.
>Mo = sum(abundance * mass)
% We also need the acceleration due to gravity (g, in m/s^2). We can't
% get accurate results by just plugging g=9.8 into the barometric
% equation because g decreases with altitude. Each number in the
% following array gives g at the corresponding altitude in the h array.
% Also assume that the pressure Po at sea level is 1 atm.
% Now calculate an array P that contains the pressures predicted
% by the barometric formula at each altitude h.
% Remember that names in Euler are case sensitive. We can use T
% for the kelvin temperature, and t for the Celsius temperature.
>T = t+273.15;
% Now use the barometric equation to calculate the total pressure,
% assuming that the average molar mass of air is the same at every
% altitude.
>title("MODEL A: Barometric Pressure Profile with Fixed Air Molar Mass");
>xlabel("Barometric Pressure (atm)");
>ylabel("altitude (m)");

% ----------------------------------------------------------------------
% TASK 2B. In MODEL A, we assumed that the average molar mass of air
% is the same at every altitude. However, each separate gas might
% have its own altitude profile, with heavier gases being more
% abundant at lower altitudes than light gases.
% On a second plot, use the barometric equation to show
% the partial pressures for each of the four gases as a function of
% altitude. Put the four curves on the same plot, in different colors.
% Call the arrays holding the partial pressures of each gas
% PN2, PO2, PAr, and PCO2.
>PN2 = (Po*abundance[1])*exp(-h*mass[1]*g/(R*T));
>PO2 = (Po*abundance[2])*exp(-h*mass[2]*g/(R*T));
>PAr = (Po*abundance[3])*exp(-h*mass[3]*g/(R*T));
>PCO2 = (Po*abundance[4])*exp(-h*mass[4]*g/(R*T));

% Then use Dalton's Law to get a fifth curve,
% the total pressure profile. Call it Ptot.
>Ptot = PN2 + PO2 + PAr + PCO2;
% Plot the curve for the total pressure with a heavier line using
% 'thickness'. For example, to plot y as a function of x with a
% heavier curve, you could type plot2d(x,y,thickness=2).
>title("MODEL B: Barometric Pressure Profile with Variable Molar Mass");
>xlabel("Pressure (atm)");
>ylabel("Altitude (m)");

% ----------------------------------------------------------------------
% TASK 2C. Calculate the variation in average molar mass with altitude
% for MODEL B. You can do this by solving the barometric equation
% for average molar mass M, using pressure Ptot (from Task 2B).
% HINT: You may run into some practical difficulties in your calculation
% of M at h = 0, because
% * division by zero is not allowed
% * ln (Ptot/Po) = 0 when h=0, so you may get molar mass = 0 at
% altitude zero.
% The average molar mass at h=0 should be Mo. You can set the second
% element of M to Mo by typing M[2] = Mo.
>h[2]=1; M = -R*T/(h*g)*ln(Ptot/Po);
Column 1 to 3:
0.02897184063066 0.02896413191 0.02895617819929
Column 4 to 6:
0.02894796718375 0.02893948906187 0.02893072545567
Column 7 to 9:
0.02892166297037 0.02891228549505 0.02890257572915
Column 10 to 12:
0.02889251507121 0.0288820868212 0.02887126708741
Column 13 to 15:
0.02882661789134 0.02878778498486 0.02875566381133
Column 16 to 18:
0.02872691182118 0.02868986042015 0.02865867484957
Column 19 to 21:
0.02859332001119 0.02852196225203 0.02845761067037

% Insert a plot the average molar mass M for model B as a function of
% altitude.
>title("MODEL B: Variation of Molar Mass with Altitude");
>xlabel("Average Molar Mass (kg/mol)");
>ylabel("Altitude (m)");

% ----------------------------------------------------------------------
% TASK 2D. Assess the accuracy of the two methods for calculating
% the absolute error in total pressure for the two methods against
% altitude, using measured pressures in the array Ptrue for the
% altitudes in h. Ptrue is in atmospheres.
>Ptrue = [1.1229908, 1.000000, 0.885790, 0.780201, 0.683079, 0.594156, 0.513204, 0.439950, 0.374090, 0.315296, 0.263226, 0.217485, 0.0949375, 0.0435330, 0.0218130, 0.0113010, 0.00455875, 0.00200189, 0.000290389, 2.35653e-5, 1.48565e-6];
>plot2d(Ptrue - Ptot, h, color=2, thickness=1);
>plot2d(Ptrue - P, h, add=1, color=3, thickness=5);
>title("Absolute Error in Model B (thinner curve) and Model A (thicker curve)");
>xlabel("Pressure (atm)");
>ylabel("Altitude (m)");

% Which of the two models is more accurate? Explain.
>"Model A is more accurate. This suggests that the gases in air don't develop separate altitude profiles.";
>"A reasonable hypothesis is that winds remix the gases that gravity would separate.";

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: