# CHEMISTRY

CHM210_A1_StratosphericChemistry

September 20, 2019

1 CHM210 Assignment 1: Stratospheric Chemistry

Welcome to the first assignment for CHM210. We are going to use a Jupyter notebook with Python to investigate some properties of stratospheric chemistry and the ozone layer. All of your an- swers will be completed within this document (changes you make are saved into your personal https://utoronto.syzygy.ca account, so don’t worry about editing this file). There are a lot of free, online resources out there to help you with Jupyter and Python. I recommend Christian Hill’s “Learning Scientific Programming with Python”. The chemistry content of this assignment should follow your textbook and lecture notes, but for those wanting additional materials, Daniel Jacob’s “Introduction to Atmospheric Chemistry” is freely available online and has several rele- vant chapters for this assignment (Ch1. “Measures of Atmospheric Composition, Ch2. ”Atmo- spheric Pressure“, and Ch10. ”Stratospheric Ozone”).

We are going to look at some real satellite data related to the structure and composition of the stratosphere in order to better understand the chemistry that takes place. The below image is data from the Atmospheric Infrared Sounder, or AIRS, instrument from NASA’s Earth Observing System (EOS) polar-orbiting platform. This data is freely available through NASA’s EarthData portal, GIOVANNI. The above data is a daily averaged vertical profile of ozone over Toronto from September 23rd, 2016. We can see where the ozone “layer” is on this graph. But why does it form there?

First things first, we need to import the Python modules we will be using for this assignment. Press the “Run” button (or ctrl-enter) on the box below to import the packages we need. If done successfully, the “You have imported the above packages!” text will display below.

1.1 Importing Python modules and getting to know Python

[2]: # lines that start with a “#” character are comments (Python doesn’t treat them␣ ↪→as code, but they are helpful for including explanations of why we have␣ ↪→written code a particular way)

# the import statements are how we add the specific packages we want to use # when we write “import YYYYYY as Y”, we are creating a shortform “Y” that we␣

↪→can use to reference the package “YYYYYY” import matplotlib.pyplot as plt import math as m import numpy as np import os

1

# below is a simple print statement that will appear as text when we run this␣ ↪→box of code

# you can easily change what text prints here by editing the words between the␣ ↪→””

print(“You have imported the above modules!”)

You have imported the above modules!

Let’s practice some simple math before we get started to make sure you’re comfortable using this Notebook. Run the code below (and feel free to change things) to see what simple arithmetic looks like in Python.

[ ]: 5+10

This is how you would write 105

[ ]: 10**5

There are several useful, pre-built, functions you can play with, like max(), min() [ ]: # max() returns the maximum value of all arguments you give it

max(1,3,6,10,5)

We can also easily write our own functions, as below: [ ]: # We are defining our own function called “mean” that will take in an “array”␣

↪→of numbers and return the mean def mean(array):

# the sum function returns the sum of all values in our array # the len function returns the number of elements in our array return sum(array) / len(array)

Now let’s use our new function [ ]: # this is just a test array (the name isn’t meaningful nor the numbers we’re␣

↪→starting with) testarray = [4,7,10,11,100] # now we call our function with the testarray as the input mean(testarray)

Now let’s look what that math package we imported allows us to do. [ ]: # the default logarithm is the natural log (ie. base e)

m.log(10)

[ ]: # this is how you call log base 10 (math.log(a,base)) m.log(10,10) # alternatively, you can use math.log10(a)

The math module contains most of the basic mathematical functions we might want, ie. sqrt(x), exp(x), sin(x), pi, e, etc. You can look them up here and test them out below:

[ ]: m.sqrt(100) #m.pi #m.exp(1)

2

1.2 Let’s load some real data

Let’s start by looking at the observed temperature profile. First, we need to download the data into our notebook. We will do so with the loadtxt function

as below. Run the code below this text box. When you run the code, the text stating “You have loaded the pressure and temperature data!” will print (and hopefully no error messages along with it).

[3]: # download .txt file containing atmospheric pressure in units of hPa Pressure = np.loadtxt(fname = “http://individual.utoronto.ca/sck/CHM210/

↪→Pressure.txt”) # download .txt file containing atmospheric temperature in units of Kelvin Temperature = np.loadtxt(fname = “http://individual.utoronto.ca/sck/CHM210/

↪→Temperature.txt”) print(“You have loaded the pressure and temperature data!”)

You have loaded the pressure and temperature data!

1.3 Learning to plot simple data

Now, let’s plot the data to see what it looks like. The code below will create a simple plot of the data we have just imported. If the below code works properly, you should see a plot with pressure on the y-axis and temperature on the x-axis.

[4]: # this makes the plot (and sets the size as 6″x8″) fig = plt.figure(figsize=(6, 8)) tp = fig.add_subplot(111)

# this tells Python what data to plot (Temperature vs Pressure) and what␣ ↪→markers to use (go– means green dashes between solid circles)

tp.plot(Temperature, Pressure, ‘go–‘) # we label the axis (with units!) tp.set_ylabel(‘Pressure (hPa)’) tp.set_xlabel(‘Temperature (K)’) # we will set the y-axis to be on a log-scale since we are plotting pressure tp.set_yscale(‘log’) # we will set the axis limits just to make things look a little nicer tp.set_ylim((1000, 1)) tp.set_xlim((210, 320))

[4]: (210, 320)

3

Now for the first question. We will walk through the first one together.

2 Q1a. Identify the height at which the tropopause occurs (you must include units with your answer!) (2 marks)

To do this, we will need to convert pressure (hPa) to altitude (km). We can do that by knowing that,

4

z = −log (

P(z) P(0)

) × (R × T)

(Mair × g) (1)

Where P(z) is the pressure at altitude z, P(0) is the pressure at the surface, R is the gas constant, T is the temperature, Mair is the molecular weight of air, and g is the acceleration due to gravity. A straightforward derivation of this can be found in Jacob. We can write this equation in Python as below:

[5]: # First, let us assign values to the variables we need to use. # When working with R, the gas constant, pay close attention to its units. R = 8.31 # m3PaK1mol1 Ma = 29 #g mol-1. g = 9.8 #m s-2

# We are going to calculate the altitude corresponding to each pressure value. # What the below code says is that for each level z, starting with the first␣

↪→pressure value (0) and ending with the last (len(Pressure)), we will␣ ↪→calculate the corresponding altitude

# for z in range(0,len(Pressure)): Altitude = [-m.log(Pressure[z]/Pressure[0])*(R*Temperature[z])/(Ma*g) for z in␣

↪→range(0,len(Pressure))]

# This will print out the contents of our array so we can make sure they make␣ ↪→sense (values should range between 0km and ~53km)

print(Altitude)

[-0.0, 0.6463918956236044, 1.3151452214061055, 2.8796906991746334, 4.011800255564502, 5.269155486326363, 6.693743805188915, 8.202566021890949, 9.130065147624338, 10.205141142150634, 11.811757392142377, 14.222340082851833, 16.3668965290799, 18.811912030935336, 22.19466832623713, 25.281169035889185, 27.193723822474034, 30.84685886963993, 33.60417303182404, 36.82317752744588, 41.8523161562467, 45.20630417570779, 48.78358518373318, 53.12719258775497]

If the above code worked, you should have printed the altitudes calculated ranging from 0km to 53km. Let’s replot our temperature vs pressure graph to show temperature vs altitude.

Fill-in the correct information in the below code (ie. replace “FILL_IN_CORRECT_NAME” with “Altitude”):

[7]: fig = plt.figure(figsize=(6, 8)) tp = fig.add_subplot(111)

tp.plot(Temperature, Pressure, ‘go–‘) tp.set_ylabel(‘Altitude (km)’) tp.set_xlabel(‘Temperature (K)’)

[7]: Text(0.5, 0, ‘Temperature (K)’)

5

Now you should be able to identify the height at which the tropopause occurs by reading it off the graph (plus or minus 5km is fine).

3 Answer to Q1a :

75km

6

4 Q1b: Is this height constant (ie. would we expect to find the tropopause at the same altitude over different parts of the world or at times of year)? Why or why not? (2 marks)

No, because increases in troposphere temperature are associated with increase in tropopause height. Factors in play include the amount of water vapour being evaporated from equatorial seas; interannual variations in tropopause height can result from both local and large-scale driv- ing forces

[double click to type your answer here – erase this text and write your answer]

5 Q2a. Plot Altitude (km) vs Number Density (molecules/cm3) (2 marks)

All axis must be labelled with correct units. To do this, you will need to remember the ideal gas law,

PV = nRT (2)

Where P, R, T are as defined above, and n is the number of moles of air contained in V volume of atmosphere. You will likely need Avogadro’s number to convert between molecules and moles.

[10]: Av = 6.023e23 # molecules mol-1 is Avogadro’s number

Write an expression to solve for the number density (Na). Check your notes or look at Jacob for a refresher if you are stuck. Look at how we wrote the expression for Altitude above to figure out how to write this piece of code. Remember to pay close attention to units!

[24]: Na = [ Av * Pressure[z]/R * Temperature[z] for z in range(0,len(Altitude))]

Once you’ve calculated number density, create your plot.

6 Answer to Q2a:

[29]: fig3 = plt.figure(figsize=(6, 8)) na = fig3.add_subplot(111)

na.plot(Na, Altitude, ‘go–‘) na.set_ylabel(‘Altitude (km)’) na.set_xlabel(‘Number Density(molecules/cm3)’)

[29]: Text(0.5, 0, ‘Number Density(molecules/cm3)’)

7

7 Q2b. Plot Altitude vs O2 (molecules/cm3) (2 marks)

All axis must be labelled with correct units.

8

8 Answer to Q2b:

[31]: # Rememeber what percentage of the atmosphere is O2 O2 = [ 0.21*Pressure[z] for z in range(0,len(Altitude))]

fig4 = plt.figure(figsize=(6, 8)) O2a = fig4.add_subplot(111) O2a.plot(O2, Altitude, ‘go–‘) O2a.set_ylabel(‘Altitude (km)’) O2a.set_xlabel(‘molecules/(cm3)’)

[31]: Text(0.5, 0, ‘molecules/cm3’)

9

9 The Chapman Mechanism

Recall from lecture the four reactions originally proposed to explain the presence of the strato- spheric ozone layer,

10

O2 + hν → O + O∗ (λ < 240nm) (R1) (3) O + O2 + M → O3 + M (R2) (4)

O3 + hν → O2 + O∗ (λ < 320nm) (R3) (5) O3 + O → 2O2 (R4) (6)

10 Q3a. Write down the rates of each of the above four reactions (2 marks)

The rate constants are k1 for R1, k2 for R2, k3 for R3, and k4 for R4. k1 and k3 have units of s−1 (photon density is already incorporated), k2 in units of cm6molecules−2s−1, and k4 in units of cm3molecules−1s−1. The rates of reactions will be written as functions of the rate constants and concentrations ([O2], [M], [O3], and [O]).

11 Answer to Q3a:

[double click to type your answer here – erase this text and write your answer]

12 Q3b. Plot the rate constants vs altitude (4 marks)

To do this, we will need to import data for k1 and k3. The data we will be using is output from the TUV calculator, a free, online tool for modelling photolysis frequencies and actinic flux.

[ ]: # download a .txt file containing photolysis rate constants for O2 -> O + O␣ ↪→from 0 – 50km~

k1 = np.loadtxt(fname = “http://individual.utoronto.ca/sck/CHM210/k1d.txt”) # download a .txt file containing photolysis rate constants for O3 -> O2 +␣

↪→O(1D) from 0 – 50km~ k3 = np.loadtxt(fname = “http://individual.utoronto.ca/sck/CHM210/k3d.txt”) print(“k1 and k3 have been loaded!”)

The expression for k2 is given by, k2 = 6 × 10−34(T/300)−2.3 You will need to write this in a form Python can interpret.

[ ]: k2 = [FILL_IN_FORMULA for z in range(0,len(Pressure))]

Similarly, k4 is given by, k4 = 8 × 10−12exp(−2060/T)

[ ]: k4 = [FILL_IN_FORMULA for z in range(0,len(Pressure))]

Now we should have everything needed to plot the four rate constants. If you have kept the names as above, the below code should plot the four rate constants vs altitude.

11

13 Answer to Q3b:

[ ]: # we are going to make two plots so k1 and k3 are on one plot and k2 and k4 are␣ ↪→on another

f, (ax1, ax2) = plt.subplots(1, 2, sharey=True) # k1 will be in red color = ‘tab:red’ ax1.set_xlabel(‘$k_1$ ($s^-1$)’, color=color) # the x-axis will be in a log-scale so we can see how many orders of magnitude␣

↪→our rate constants span ax1.set_xscale(‘log’) ax1.set_ylabel(‘Altitude (km)’) ax1.plot(k1, Altitude, color=color) ax1.tick_params(axis=’x’, labelcolor=color) ax1.set_xlim((1e-16, None))

ax3 = ax1.twiny() # k3 will be in red color = ‘tab:blue’ ax3.set_xlabel(‘$k_3$ ($s^-1$)’, color=color) ax3.set_xscale(‘log’) ax3.plot(k3, Altitude, color=color) ax3.tick_params(axis=’x’, labelcolor=color)

# k2 will be in green color = ‘tab:green’ ax2.set_xlabel(‘$k_2$ ($cm^{6}molecules^{-2}s^{-1}$)’, color=color) ax2.plot(k2, Altitude, color=color) ax2.tick_params(axis=’x’, labelcolor=color) ax2.set_xlim((0, None))

ax4 = ax2.twiny() # k4 will be in cyan color = ‘tab:cyan’ ax4.set_xlabel(‘$k_4$ ($cm^{3}molecules^{-1}s^{-1}$)’, color=color) ax4.plot(k4, Altitude, color=color) ax4.tick_params(axis=’x’, labelcolor=color) ax4.set_xlim((0, None))

12

14 Q3c. What is the steady-state assumption and when might it apply to chemical species in the atmosphere (2 marks)

15 Answer to Q3c:

[double click to type your answer here – erase this text and write your answer]

16 Q3d. Write an expression for the steady-state concentration of O3 as a function of the rate constants and [O2] only. (4 marks)

17 Answer to Q3d:

[double click to type your answer here – erase this text and write your answer]

18 Q3e. Plot the the steady-state concentration of O3 vs altitude (2 marks)

19 Answer to Q3e:

[ ]: # You will need to convert your expression for the steady-state concentration␣ ↪→of ozone into a form Python can recognize

O3 = [FILL_IN_FORMULA for z in range(0,len(Pressure))]

fig5 = plt.figure(figsize=(6, 8)) O3a = fig5.add_subplot(111)

O3a.plot(O3, Altitude, ‘go–‘) O3a.set_ylabel(‘Altitude (km)’) O3a.set_xlabel(‘[$O_3$] (molecules/$cm^3$)’)

13

20 Q3f. By referencing the above plots you have made, explain why the ozone layer forms where it does. (3 marks)

21 Answer to Q3f.

[double click to type your answer here – erase this text and write your answer]

22 Q4a. Plot our steady-state [O3] vs altitude on the same graph as the satellite observations above in units of ppbv. (3 marks)

We need to convert our data from molecules/cm3 to ppbv. We will also have to download the ozone observations in a form we can easily use:

[ ]: # download .txt containing [O3] observations from 0 – 50km~ O3obs = np.loadtxt(fname = “http://individual.utoronto.ca/sck/CHM210/O3.txt”)

Now convert from molecules/cm3 to ppb [ ]: # R = 82.057338 # cm3 atm K=1 mol=1

O3ppbv = [FILL_IN_FORMULA for z in range(0,len(Pressure))]

23 Answer to Q4a:

[ ]: fig5 = plt.figure(figsize=(6, 8)) O3p = fig5.add_subplot(111)

# this will plot both curves on the same graph # our calculate O3 will be in green and the observations in red O3p.plot(O3ppbv, Pressure, ‘go–‘, label=’Modelled Ozone’) O3p.plot(O3obs, Pressure, ‘r+-‘, label=’Observed Ozone’) O3p.legend(loc=’upper right’) O3p.set_ylabel(‘Pressure (hPa)’) O3p.set_xlabel(‘[$O_3$] (ppbv)’) O3p.set_yscale(‘log’) O3p.set_ylim((1000, 1))

14

24 Q4b. Determine the overhead ozone column depth in Dobson units for both our calculated O3 and the observed O3 (DU) (4 marks)

To calculate the ozone column depth, we need to integrate the total ozone vs altitude. This might be easier to think about if we flip out axis, ie.

[ ]: fig6 = plt.figure(figsize=(6, 8)) O3p = fig6.add_subplot(111) O3p.plot(Altitude, O3, ‘r+-‘, label=’Observed Ozone’) O3p.set_xlabel(‘Altitude (km)’) O3p.set_ylabel(‘[$O_3$] (molecules/$cm^3$)’)

We want to find the area under the curve above. We can do that with the numpy.trapz func- tion, which uses the trapezoidal rule for numerical integration. We’ll need to make sure units are appropriate (ie. altitude should be in cm if concentration is in molecules/cm3).

[ ]: Altitudecm = [Altitude[z]*1e5 for z in range(0,len(Pressure))] np.trapz(O3, x=Altitudecm)

25 Answer to Q4b.

[double click to type your answer here – erase this text and write your answer]

26 Q4c. Why might modelled ozone be so much greater than the ob- served ozone concentrations? (4 marks)

27 Answer to Q4c.

[double click to type your answer here – erase this text and write your answer]

When you have finished your assignment, you need to save it as a PDF and submit that PDF on Quercus. In your Jupyter Notebook, goto File –> Download ad –> PDF via LaTeX (.pdf)

If you have technical questions related to Python or Jupyter, contact Sarah Kavassalis (sarah.kavassalis@mail.utoronto.ca).

15

- CHM210 Assignment 1: Stratospheric Chemistry
- Importing Python modules and getting to know Python
- Let’s load some real data
- Learning to plot simple data

- Q1a. Identify the height at which the tropopause occurs (you must include units with your answer!) (2 marks)
- Answer to Q1a :
- Q1b: Is this height constant (ie. would we expect to find the tropopause at the same altitude over different parts of the world or at times of year)? Why or why not? (2 marks)
- Q2a. Plot Altitude (km) vs Number Density (molecules/cm^3) (2 marks)
- Answer to Q2a:
- Q2b. Plot Altitude vs O_2 (molecules/cm^3) (2 marks)
- Answer to Q2b:
- The Chapman Mechanism
- Q3a. Write down the rates of each of the above four reactions (2 marks)
- Answer to Q3a:
- Q3b. Plot the rate constants vs altitude (4 marks)
- Answer to Q3b:
- Q3c. What is the steady-state assumption and when might it apply to chemical species in the atmosphere (2 marks)
- Answer to Q3c:
- Q3d. Write an expression for the steady-state concentration of O_3 as a function of the rate constants and [O_2] only. (4 marks)
- Answer to Q3d:
- Q3e. Plot the the steady-state concentration of O_3 vs altitude (2 marks)
- Answer to Q3e:
- Q3f. By referencing the above plots you have made, explain why the ozone layer forms where it does. (3 marks)
- Answer to Q3f.
- Q4a. Plot our steady-state [O_3] vs altitude on the same graph as the satellite observations above in units of ppbv. (3 marks)
- Answer to Q4a:
- Q4b. Determine the overhead ozone column depth in Dobson units for both our calculated O_3 and the observed O_3 (DU) (4 marks)
- Answer to Q4b.
- Q4c. Why might modelled ozone be so much greater than the observed ozone concentrations? (4 marks)
- Answer to Q4c.