MODELLING PHOTOCHEMICAL AIR POLLUTANTS FROM INDUSTRIAL ANDBIOGENIC EMISSIONS IN THE TERRACE - KITIMAT VALLEY:A CONSTRAINED COASTAL AIRSHED WITH COMPLEX TERRAINbyBenjamin Ralph WeinsteinB.Sc, The University of British Columbia, 2003A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Atmospheric Science)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2015© Benjamin Ralph Weinstein, 2015AbstractThis study investigates the worst-case ambient concentrations of O3 and its precursors (NOXand VOCs) that may occur from the construction of large industrial facilities in the Terrace-Kitimat valley airshed. This research is important as the Terrace-Kitimat valley naturallyemits high levels of biogenic VOCs in the summer and many of the proposed facilities will, ifconstructed, emit high quantities of NOX . To date, literature concerning O3 production formindustrial development in coastal airsheds with complex terrain is sparse.The Comprehensive Air Quality Model with Extensions (CAMx) was used as the photochem-ical model for this research. Spring and summer periods were selected from 2010. Controland Test Case emission inventories were developed, the former for model evaluation and thelatter to assess pollutant change. Model evaluation showed that CAMx was able to emulatedaytime O3 peaks in an adjacent valley for both periods though overnight titration by NOwas less adequately replicated. Sensitivity tests revealed that this was due in part to inade-quate Control Case emissions quantification; results improved with the addition of small scalearea-based NOX emissions to account for missing sources in the original emissions inventory.Results from the spring period suggest that increased industrial emissions, in general, wouldnot contribute to valley-wide O3 increases greater than 5 ppb, as biogenic VOC emissionsare minimal throughout the airshed during this season. On the other hand, results from thesummer period suggest that increased industrial emissions would, at times, contribute to agreater than 55% increase in O3 concentrations, particularly downwind of Kitimat on dayswith high temperatures, low planetary boundary layer heights, differential heating of the landand ocean surface temperatures and consecutive days of horizontally recirculating wind.This research also used the modelled O3 - reactive nitrogen ratio during hours conducive tophotochemistry to determine the O3 sensitivity of the Terrace-Kitimat valley airshed. Theiiairshed is currently sensitive to NOX emissions however the full construction of all proposedindustrial projects would likely change the O3 sensitivity of a large portion of the valley to besensitive to emissions of VOCs, especially in and around Kitimat.iiiPrefaceThis thesis contains the original research and analysis undertaken by the author BenjaminRalph Weinstein, under the guidance of supervisors Douw Steyn and Peter Jackson.This thesis contains many figures, some of which are maps. All maps were created by the authorusing the computer software QGIS, freely available for download from http://www.qgis.org/en/site/. Map base layers were obtained from: online sources (http://geogratis.cgdi.gc.ca/), Dave Amirault and Blair Ells at the B.C. Ministry of Forests, Lands and NaturalResource Operations, and Morgan Hite at http://hesperus-wild.org/.For the purpose of this research WRF output was obtained from David Suita, a PhD studentworking at the Geophysical Disaster Computational Fluid Dynamics Centre at UBC.Results from Chapters 3 and 4 of this thesis were presented at the International TechnicalMeeting on Air Pollution Modelling and its Application in Montpellier, France, in March 2015.The title of this presentation was: Modelling photochemical air pollutants from industrialemissions in a constrained coastal valley with complex terrain. The extended abstract from thispresentation is forthcoming in: Air Pollution Modeling and its Application XXIV, publishedby Springer.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvList of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Ozone and Ozone Chemistry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 Production of Ozone in Polluted Air . . . . . . . . . . . . . . . . . . . . . . 21.1.2 Canadian Standards and Objectives for Ground Level Ozone . . . . . . . 41.2 Proposed Industrial Development in the Terrace - Kitimat Valley . . . . . . . . . 61.2.1 Liquefied Natural Gas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.3 Photochemical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.3.1 Comprehensive Air Quality Model with Extensions . . . . . . . . . . . . . 10v1.4 Research Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.5 Proposed Research Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.1 Period Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2 Developing Meteorological Inputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.2.1 Weather Research and Forecast Model . . . . . . . . . . . . . . . . . . . . . 192.2.2 Weather Research and Forecast Model Evaluation . . . . . . . . . . . . . . 232.2.3 Intermediate Programs: MCIP and WRFCAMx . . . . . . . . . . . . . . . 262.3 Developing Emissions Inputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.3.1 Anthropogenic Emissions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.3.2 Biogenic Emissions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.3.3 Other Emissions Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . 482.4 Developing Boundary and Initial Conditions Inputs . . . . . . . . . . . . . . . . . 502.5 Developing Other Model Inputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 502.6 CAMx Model Settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 512.6.1 Advection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 512.6.2 Chemical Mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 512.6.3 Aerosol Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 512.6.4 Plume in Grid Option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 522.6.5 Wet and Dry Deposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533 Results and Discussion I: Answering Question 1 . . . . . . . . . . . . . . . . . . 543.1 Model Fitness Evaluation Description . . . . . . . . . . . . . . . . . . . . . . . . . 553.2 Evaluation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.2.1 Ozone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.2.2 Fitness of Control Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.3 CAMx Sensitivity to Mobile and Point Source Emissions Perturbations . . . . . 613.3.1 Sensitivity Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61vi3.3.2 Ozone Sensitivity Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 633.4 Question 1 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . 674 Results and Discussion II: Answering Question 2 . . . . . . . . . . . . . . . . . 714.1 Springtime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.1.1 Time Series: Communities in the Terrace - Kitimat Valley Airshed . . . 724.1.2 Locations of Minima and Maxima . . . . . . . . . . . . . . . . . . . . . . . 734.1.3 Ingredients for Reduced Ozone . . . . . . . . . . . . . . . . . . . . . . . . . 784.2 Summertime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.2.1 Time Series: Hotspots in the Terrace - Kitimat Valley Airshed . . . . . . 794.2.2 Locations of Minima and Maxima . . . . . . . . . . . . . . . . . . . . . . . 804.2.3 Ingredients for Elevated Ozone . . . . . . . . . . . . . . . . . . . . . . . . . 854.2.4 Effect of a Wildfire on Ozone Production . . . . . . . . . . . . . . . . . . . 965 Results and Discussion III: Answering Question 3 . . . . . . . . . . . . . . . . .1045.1 Method for Determining NOX -VOC Sensitivity . . . . . . . . . . . . . . . . . . . . 1055.2 Sensitivity Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1076 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .1116.1 Findings: Research Questions Revisited . . . . . . . . . . . . . . . . . . . . . . . . 1116.2 Future Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .116Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .126A Random Forest Regression Method for Modelling Ozone in Kitimat . . . .126A.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126A.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127A.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128B Statistical Metrics Used for Model Evaluation . . . . . . . . . . . . . . . . . . .131viiC WRF Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .133C.1 Quantitative Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133C.2 Qualitative Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144C.3 WRF Evaluation Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . 151D Control Case Model Evaluation for NO and NO2 . . . . . . . . . . . . . . . . .152D.1 Nitric Oxide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152D.2 Nitrogen Dioxide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154D.3 Vertical layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155E Sensitivity Test Results for NO and NO2 . . . . . . . . . . . . . . . . . . . . . . .158E.1 Nitric Oxide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158E.2 Nitrogen Dioxide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161E.3 Scatterplots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163viiiList of Tables1.1 Canadian Ambient Air Quality Standards for O3 . . . . . . . . . . . . . . . . . . . 51.2 Status of proposed projects in Kitimat . . . . . . . . . . . . . . . . . . . . . . . . . 72.1 Dates for the spring and summer model periods (2010). . . . . . . . . . . . . . . . 192.2 Properties of the WRF model domains . . . . . . . . . . . . . . . . . . . . . . . . . 202.3 Physics options used in WRF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.4 Properties of the meteorological monitoring stations used for the WRF evaluation. 232.5 Parameters measured at meteorological monitoring stations in the Terrace -Kitimat valley airshed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.6 CAMx input models sorted by translator program . . . . . . . . . . . . . . . . . . 262.7 Properties of the CAMx model domains . . . . . . . . . . . . . . . . . . . . . . . . 272.8 Biogenic volatile organic compound species output by the Model of Gases andAerosols from Nature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.1 Location of Smithers St. Joseph air quality monitoring station and correspond-ing x and y-grid cells in the 4.0 km domains. . . . . . . . . . . . . . . . . . . . . . 553.2 Dates removed from model evaluation . . . . . . . . . . . . . . . . . . . . . . . . . 573.3 O3 model evaluation statistics for the Control case - (all hours) . . . . . . . . . . 593.4 O3 model evaluation statistics for the Control case - (10:00 to 21:00) . . . . . . . 593.5 Sensitivity test names and totalled NOX emissions in Smithers . . . . . . . . . . 623.6 O3 model evaluation statistics for the sensitivity tests - (all hours) . . . . . . . . 653.7 O3 model evaluation statistics for the sensitivity tests - (10:00 to 21:00) . . . . . 654.1 Dates meeting screening criteria for elevated O3 in the Terrace - Kitimat valley 95ix5.1 Hours used for O3/NOy analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107A.1 Model evaluation statistics for spring and summer periods at Smithers: multiplelinear regression model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128A.2 Model evaluation statistics for spring and summer periods at Smithers: neuralnetwork regression model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128A.3 Model evaluation statistics for spring and summer periods at Smithers: randomforest regression model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128A.4 Model evaluation statistics for spring and summer periods at Kitimat usingmodel coefficients developed from Smithers: random forest regression model . . 129A.5 Dates for the spring and summer model periods (2010). . . . . . . . . . . . . . . . 130C.1 Modelled - observed temperature statistics from the WRF evaluation . . . . . . 134C.2 Modelled - observed RH statistics from the WRF evaluation . . . . . . . . . . . . 136C.3 Modelled - observed pressure statistics from the WRF evaluation . . . . . . . . . 137C.4 Modelled - observed wind direction statistics from the WRF evaluation . . . . . 140C.5 Modelled - observed wind speed statistics from the WRF evaluation . . . . . . . 141D.1 NO model evaluation statistics for the Control case - (all hours) . . . . . . . . . . 153D.2 NO2 model evaluation statistics for the Control case - (all hours) . . . . . . . . . 154E.1 NO model evaluation statistics for the sensitivity tests - (all hours) . . . . . . . . 159E.2 NO2 model evaluation statistics for the sensitivity tests - (all hours) . . . . . . . 161xList of Figures1.1 Canadian air management threshold values and actions. . . . . . . . . . . . . . . 61.2 TKVA overview with 250 m contours. . . . . . . . . . . . . . . . . . . . . . . . . . 81.3 A broad schematic of the Eulerian photochemical grid modelling system. . . . . 111.4 The overall modelling framework and input categories in CAMx. . . . . . . . . . 122.1 Maximum daily O3 mixing ratios at Smithers. . . . . . . . . . . . . . . . . . . . . 182.2 All WRF domains. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.3 Outline and individual cell structure of the WRF d04. . . . . . . . . . . . . . . . 222.4 Evaluation locations for WRF d04 output. . . . . . . . . . . . . . . . . . . . . . . . 242.5 Outlines of WRF d04 and CAMx grd03. . . . . . . . . . . . . . . . . . . . . . . . . 282.6 The vertical levels in WRF mapped to the vertical levels in CAMx. . . . . . . . 292.7 Emission locations in the TKVA for the Control case. . . . . . . . . . . . . . . . . 312.8 Daily CAC emissions in the TKVA for the Control case. . . . . . . . . . . . . . . 322.9 BV and TKV overview and emission locations in the BV for the Control case. . 342.10 Daily CAC emissions in the BV for the Control case. . . . . . . . . . . . . . . . . 352.11 Emission locations in the TKVA for the Test case. . . . . . . . . . . . . . . . . . . 372.12 Daily CAC emissions in the TKVA for the Test case. . . . . . . . . . . . . . . . . 382.13 Biogeoclimatic zones and landcover types of the TKV. . . . . . . . . . . . . . . . 402.14 Daily bVOC emissions in the TKV airshed. . . . . . . . . . . . . . . . . . . . . . . 422.15 Isoprene emissions in the TKV on August 15th, 2010 at 16:00 PST. . . . . . . . . 442.16 Biogeoclimatic zones and landcover types of the BV. . . . . . . . . . . . . . . . . 462.17 Daily bVOC emissions in the BV airshed. . . . . . . . . . . . . . . . . . . . . . . . 47xi2.18 Isoprene emissions in the BV and TKV. . . . . . . . . . . . . . . . . . . . . . . . . 493.1 Observed and modelled O3 time series at Smithers for the spring and summermodel periods, Control Case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583.2 Observed vs. modelled O3 scatterplot at Smithers for spring and summer modelperiods, Control Case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.3 Hourly distribution of daily point and mobile emissions. . . . . . . . . . . . . . . 633.4 O3 sensitivity test time series for spring and summer model periods. . . . . . . . 643.5 Observed vs. modelled O3 scatterplot at Smithers for spring and summer modelperiods, 20xp sensitivity case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663.6 Updated daily NOX emissions in the TKV for the Test case. . . . . . . . . . . . . 693.7 Updated emission locations in the TKVA for the Test case. . . . . . . . . . . . . 704.1 PBL depth, O3 and temperature time series for the spring model period atTerrace, Kitimat and Hartley Bay. . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.2 CAMx spatial output for the TKVA on April 17th, 2010 at 00:00 PST showingTest Case modelled O3 mixing ratios, and percent change from the Control case. 754.3 CAMx spatial output for the TKVA on May 4th, 2010 at 16:00 PST showingTest Case modelled O3 mixing ratios, and percent change from the Control case. 774.4 PBL depth time series for the spring model period in the TKVA. . . . . . . . . . 784.5 PBL depth, O3 and temperature time series for the summer model period atLakelse Lake, Kitimat and Miskatla Inlet. . . . . . . . . . . . . . . . . . . . . . . . 814.6 CAMx spatial output for the TKVA on August 4thth, 2010 at 15:00 PST show-ing Test Case modelled O3 mixing ratios, and percent change from the Controlcase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.7 CAMx spatial output for the TKVA on August 15th, 2010 at 16:00 PST showingTest Case modelled O3 mixing ratios, and percent change from the Control case. 844.8 Maximum daily temperature vs. bVOC emissions scatterplots in the Terrace -Kitimat valley airshed for spring and summer model periods. . . . . . . . . . . . 86xii4.9 PBL depth time series for the summer model period in the TKVA. . . . . . . . . 874.10 Difference between Kitimat and Douglas Channel 2m temperature. . . . . . . . . 884.11 Wind streamlines and wind barb cross sections in the TKV on August 14th. . . 894.12 Wind streamlines and wind barb cross sections in the TKV on August 15th . . . 904.13 CAMx spatial O3 output for the TKVA on August 14th, 2010. . . . . . . . . . . . 924.14 CAMx spatial O3 output for the TKVA on August 15th, 2010. . . . . . . . . . . . 934.15 Haze in the BV from a nearby forest fire on August 14th and August 15th. . . . 974.16 Forest fire locations within the 12 km model domain on August 14th 2010. . . . 984.17 Spatial NO fields generated by the MOZART for use as initial conditions toCAMx for July 29th, August 12th, August 13th and August 14th. . . . . . . . . . 994.18 Spatial NO2 fields generated by the MOZART for use as initial conditions toCAMx for July 29th, August 12th, August 13th and August 14th. . . . . . . . . . 1004.19 Spatial O3 fields generated by the MOZART for use as initial conditions toCAMx for July 29th, August 12th, August 13th and August 14th. . . . . . . . . . 1004.20 Lakelse Lake O3 time series with varying ICs . . . . . . . . . . . . . . . . . . . . . 1014.21 CAMx spatial O3 output for the TKVA on August 15th, 2010 at 16:00 PST forTest Case model initialization on July 29th and August 14th. . . . . . . . . . . . . 1035.1 North - south transect of the TKV. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1065.2 Median, mean, 5th and 95th percentile O3/NOy ratios along the TKV transect. 1085.3 Median O3/NOy ratios in the TKVA for the Control and Test cases. . . . . . . . 110A.1 Modelled O3 time series at Kitimat using RF model coefficients developed inSmithers for spring and summer model periods. . . . . . . . . . . . . . . . . . . . 130C.1 Observed and modelled temperature time series at Terrace, CYXT and Kitimatfor spring and summer model periods. . . . . . . . . . . . . . . . . . . . . . . . . . 134C.2 Observed and modelled RH time series at Terrace and CYXT for spring andsummer model periods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136xiiiC.3 Observed and modelled pressure time series at CYXT for spring and summermodel periods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138C.4 Observed and modelled wind direction time series at Terrace, CYXT and Kiti-mat for spring and summer model periods. . . . . . . . . . . . . . . . . . . . . . . 139C.5 Observed and modelled wind speed time series at Terrace, CYXT and Kitimatfor spring and summer model periods. . . . . . . . . . . . . . . . . . . . . . . . . . 143C.6 Wind streamlines in the TKV on July 31st, 2010 at 00:00, 04:00, 08:00 and12:00 PST. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146C.7 Wind streamlines in the TKV on July 31st, 2010 at 16:00, 20:00 and 00:00(August 1st) PST. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147C.8 Temperature evolution for the morning of July 31st, 2010. . . . . . . . . . . . . . 149C.9 Temperature evolution for the afternoon of July 31st, 2010. . . . . . . . . . . . . 150D.1 Observed and modelled NO time series at Smithers for the spring and summermodel periods, Control Case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153D.2 Observed and modelled NO2 time series at Smithers for the spring and summermodel periods, Control Case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155D.3 PBL depth, wind speed and vertical profiles of NO2 and NO at Smithers fromMay 6th - 12th, 2010. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157E.1 NO sensitivity test time series for spring and summer model periods. . . . . . . . 160E.2 NO2 sensitivity test time series for spring and summer model periods. . . . . . . 162E.3 Spring: modelled vs. observed O3, NO and NO2. . . . . . . . . . . . . . . . . . . . 164E.4 Summer: modelled vs. observed O3, NO and NO2. . . . . . . . . . . . . . . . . . . 165xivList of Symbols↵ one of two isomers of pinene wavelengthµ microxvList of AcronymsAGL above ground levelAQM Air Quality ModelAQMS Air Quality Management SystemBC British ColumbiaBV Bulkley ValleyBVLD Bulkley Valley - Lakes DistrictbVOC biogenic VOCsCAAQS Canadian Ambient Air Quality StandardsCAC criteria air contaminantsCAMS Community Modelling and Analysis SystemCAMx Comprehensive Air Quality Model with ExtensionsCB05 Carbon Bond (2005 update)CF Coarse-FineCH4 methaneCMAQ Community Multiscale Air Quality Modelling SystemCMU Carnegie Mellon UniversityCO carbon monoxidexviCREATE-AAP Collaborative Research and Training Experience - Atmospheric AerosolProgramGREASED Greatly Reduced Execution and Simplified DynamicsICs Initial conditionsLFV Lower Fraser ValleyLNG liquefied natural gasMAE mean absolute errorMB mean biasMCIP Meteorology-Chemistry Interface ProcessorMEGAN Model of Emissions of Gases and Aerosols from NatureMLR multiple linear regressionMOE Ministry of EnvironmentMOZART Model for Ozone and Related Chemical TracersNAAQO National Ambient Air Quality ObjectivesNAM North American Mesoscale Forecast SystemNCAR National Centre for Atmospheric ResearchNCEP National Centre for Environmental PredictionNCL NCAR Command LanguageNMB normalized mean biasNMHC non-methane hydrocarbonsxviiNN neural networksNO nitric oxideNO2 nitrogen dioxideNOX oxides of nitrogenNOAA National Oceanic and Atmospheric AdministrationNOy total reactive nitrogenNPRI National Pollutant Release InventoryNSERC Natural Sciences and Engineering Research Council of CanadaO3 OzoneOH hydroxyl radicalPAN peroxyacetyl nitratePBLH planetary boundary layer heightPiG Plume-in-GridPM particulate matterPM2.5 particulate matter with aerodynamic diameter less than 2.5 micrometersppb parts per billionPPM Piecewise Parabolic Methodppm parts per millionPST Pacific Standard Timer correlation coefficientRF random forestsxviiiRH relative humidityRMSE root mean squared errorRTA Rio Tinto AlcanSMOKE Sparse Matrix Operator Kernel Emissions Modelling SystemSO2 sulphur dioxideTKV Terrace - Kitimat ValleyTKVA Terrace - Kitimat Valley AirshedVOCs volatile organic compoundsWRF Weather Research and Forecast ModelxixAcknowledgementsI would like to begin this section by recognizing Douw Steyn. Douw’s commitment to hisstudents far exceeds what could be considered the normal requirements of a supervisor. Douw’smentorship, friendship, encouraging nature and overall support enabled me to complete thiswork, and for this I am truly grateful. I would also like to thank my other supervisor, PeterJackson, for his support along the way and Ian McKendry who is the third member of mygraduate committee.I gratefully acknowledge the Natural Sciences and Engineering Research Council of Canada(NSERC) Collaborative Research and Training Experience - Atmospheric Aerosol Program(CREATE-AAP) for funding large portions of my MSc. research and travel, and for an excel-lently structured program that brings all kinds of aerosol scientists together. The remainderof my research was funded through NSERC Discovery grants to Douw Steyn for which I amgrateful.I would also like to thank colleagues at the B.C. Ministry of Environment, specifically IanSharpe, Ed Hoffman and Greg Tamblyn for allowing me to take a two-year educational leave,as well as Arvind Saraswat and Robyn Roome for supporting me after I left. I would like torecognize the Pacific Leaders program for providing me with funding for tuition and books inmy second year.I would like to recognize Jeff Lundgren from RWDI Inc. who provided me with valuabletechnical support and insights as I struggled against ’big data’, as well as Morgan Hite fromHesperous Wild who taught me GIS and helped with the analysis of mobile source emissions.Special thanks to my lab mates, Annie S. and Vlad P. for the entertaining dialogue andcompany while I lived in Vancouver, and to Nadya M. for the discussions in France.xxDedicationMost importantly, I would like to thank my wife Jessica Hall. Jess, without you this would stillbe a far off dream. You have supported me throughout this journey, allowing our family to beuprooted and relocated to Vancouver. You managed our affairs while I studied and allowed meto become an additional dependent, particularly when we lived in Vancouver. Your constantsupport enabled me to succeed. For that I will always be grateful.xxiChapter 1Introduction1.1 Ozone and Ozone ChemistryOzone (O3), a colourless molecule containing three oxygen atoms, plays two distinct roles inEarth’s atmosphere. In the stratosphere (approx 15 to 50 km above ground level (AGL))O3 is a valued greenhouse gas, protecting Earth’s surface from biologically-damaging UV-Bradiation (Vallero, 2014). In the troposphere however (surface to 15 km), exposure to O3can be harmful to human health and vegetation (Krupa, 1997). O3 is most plentiful in thelower part of the stratosphere commonly referred to as the O3 layer (approximately 15-28 km)(Finlayson-Pitts and Pitts Jr, 2000). International treaties such as the Montreal Protocol andthe Vienna Convention have been developed in order to protect and restore O3 in this part ofthe atmosphere (see: Ozone Secretariat, 2015), while at Earth’s surface government standardsexist to protect people and plants from exposure to O3 in high concentrations.O3 is not emitted directly into the troposphere; it is either transported from the stratosphereor formed through the chemical reactions of its precursors in the presence of sunlight. Forthis reason O3 is considered a secondary pollutant as opposed to a primary pollutant (thosepollutants emitted directly into the atmosphere). Tropospheric O3 chemistry is complex,1nonlinear and involves a multitude of competing reactions. The majority of tropospheric O3is produced through chemical reactions (3000 - 4600 Tg⋅O3⋅ yr−1), while a smaller amountis stratospheric in origin (400 - 1100 Tg⋅O3⋅ yr−1) (Jacob, 1999). Tropospheric O3 is largelyproduced from the photolysis of nitrogen dioxide (NO2):NO2 + hv O2 NO + O3 (<420 nm) {1.1}NO2 is often combined with nitric oxide (NO) and referred to as oxides of nitrogen (NOX).These gases have a variety of tropospheric origins, the largest of which involve high temperaturecombustion (fossil fuels, biomass burning, etc.) (Jacob, 1999). In many respects, the questionof O3 formation is less about what is making O3 and more about what is making NO2.O3 photochemistry exhibits a cycling nature with NO and NO2 based on the following returnreaction:NO + O3 NO2 + O2 {1.2}This cycle leads to a photostationary state where the O3 concentration can be determinedfrom the following expression: [O3] = J3 ⋅ [NO2]k2 ⋅ [NO] (1.1)where J 3 is the rate of NO2 photolysis and k2 is the rate coefficient for reaction 1.2 (Jenkinand Clemitshavm, 2002). At night when there is insufficient solar radiation this balance shiftstowards reaction 1.2, and O3 concentrations decrease as a result (Sillman, 2012).1.1.1 Production of Ozone in Polluted AirIn order to produce tropospheric O3 in elevated concentrations, the presence of other gasesis necessary. Interestingly, reactions that eventually produce O3 all start with O3 moleculesthemselves. O3 is a key precursor to the hydroxyl radical (OH), an important molecule that2controls the oxidizing power of the troposphere through the reactions:O3 + hv O(1D) + O2 (<325 nm) {1.3}O(1D) + H2O 2OH {1.4}(Chin et al., 1994)Once liberated, OH attacks many otherwise inert gasses such as carbon monoxide (CO) andmethane (CH4). CO is a product of incomplete combustion (typically of fossil fuels) (Fenger,2002), while CH4, a powerful greenhouse gas, has a variety of anthropogenic (natural gasleakage, coal mining and petroleum industries, landfills and the raising of livestock) and naturalsources (wetlands are the largest) (U.S. Environmental Protection Agency, 2015a). In thepresence of NOX , oxidation of CO and CH4 can lead to modest production of additional O3(Jacob, 1999).In the presence of NOX and non-methane hydrocarbons (NMHC), OH chemistry can lead tobroad O3 production and photochemical smog (Crutzen, 1995). Basic reactions involved withO3 air pollution are presented below, while those concerning CO and CH4 can be found in theliterature (succinctly summarized in (Crutzen, 1995) and (Crutzen, 1979)).Following the production of OH via reactions 1.3 and 1.4, where where NOX and NMHCabound (NMHC are a subset of volatile organic compounds (VOCs)), additional NO2 forma-tion takes place as part of a five-sequence chain.RH + OH R + H2O {1.5}R + O2 + M RO2 + M {1.6}RO2 + NO RO + NO2 {1.7}RO + O2 R′CHO + HO2 {1.8}HO2 + NO OH + NO2 {1.9}3Where the term RH denotes an organic compound, RO2 denotes the subsequent peroxy radical,RO is a monoxy radical and R’CHO denotes a radical aldehyde. Considering that each NO2molecule is subsequently photolized into NO and O3 via reaction 1.1, the net reaction for thissequence is:RH + 4O2 R′CHO + 2O3 + H2O {1.10}O3 production is complex and nonlinear. In areas with low NOX , O3 production varies lin-early with increasing NO concentrations, while in areas with high NOX , O3 production varieslinearly with increasing hydrocarbon concentrations but inversely with NO2 concentrations(Jacob, 1999). Thus the emission of one precursor does not guarantee elevated ozone; at timesincreasing NOX emissions leads to reduced O3 concentrations.The cycling in reactions 1.5 through 1.9 ceases when the supply of HOX (OH + HO2) runsout. This can happen in one of two ways:HO2 + HO2 H2O2 + O2 {1.11}NO2 + OH + M HNO3 + M {1.12}Reaction 1.11 takes place when NOX concentrations are low, while 1.12 takes place when NOXconcentrations are very high (Jacob, 1999).1.1.2 Canadian Standards and Objectives for Ground Level OzoneAcute exposure to elevated O3 has been linked with negative health outcomes such as morbidityand mortality (Tiwary and Colls, 2010). In order to protect the public from exposure toO3 (along with other air pollutants such as particulate matter with aerodynamic diameterless than 2.5 micrometers (PM2.5)), the Canadian Council of Ministers of the Environment(CCME) has developed the Air Quality Management System (AQMS). This collaborativeapproach for improving air quality in Canada contains many elements, including: the CanadianAmbient Air Quality Standards (CAAQS), the concept of airsheds and air zones, industrial4emission standards, mobile source emission standards and monitoring/reporting requirements(Canadian Council of Ministers of the Environment, 2012b). Both the federal and provincialgovernments have roles in the implementation of the AQMS, outlined by the Canadian Councilof Ministers of the Environment (2012c).As part of the AQMS the federal government has developed the CAAQS for O3 and PM2.5.The O3 standard can be found in Table 1.1. In addition to the CAAQS, the British Columbia(BC) government also has a 1-hour advisory threshold of 82 parts per billion (ppb) which isbased on an older set of national objectives, the National Ambient Air Quality Objectives(NAAQO) (British Columbia Ministry of Environment, 2014).Table 1.1: Canadian Ambient Air Quality Standards for O3Pollutant Averaging Time Standard Metric2015 2020O3 8-hour 63 ppb 62 ppb 3-year average of the annual4th-highest daily maximum 8-hour average concentrations.Within the context of the AQMS the above standards is not the starting point for air qualitymanagement purposes; a tiered approach is outlined by the Canadian Council of Ministers ofthe Environment (2012b) whereby a series of management actions are meant to occur as airquality begins to deteriorate (based on monitoring results).These tiered levels, along with their respective ambient values are presented in figure 1.1. Themetric associated with these threshold values is the same as the O3 metric for the CAAQS(Canadian Council of Ministers of the Environment, 2012a), and is shown in table 1.1. Theintent of this process is to keep clean areas clean and provide motivation for air quality man-agement well before the CAAQS is exceeded.5Figure 1.1: Canadian air management threshold values and actions (Canadian Council ofMinisters of the Environment, 2012b).1.2 Proposed Industrial Development in the Terrace - KitimatValleyOne area which may require some form of air quality management related to the transitioningfrom the lowest category (green) to some other coloured category (unknown at this time) is theTerrace - Kitimat Valley (TKV), located at the head of the Douglas channel in central-west BC.The TKV stretches 55 km from Kitimat in the south to Lakelse Lake in the north and includescomplex coastal terrain with a narrow valley bottom at sea level and mountains reaching 1500m in elevation. Within the TKV is the District of Kitimat, pop. 8,300 (Statistics Canada,2014a). To the southeast of Kitimat is Kitamaat Village, a small Haisla community. Justoutside the northern edge of the TKV is the city of Terrace which lies along the Skeena river.Just east of Terrace is Thornhill, an unincorporated community. Together these communitieshave a population of approximately 11,500 (Statistics Canada, 2015a).The TKV is the primary valley in the Terrace - Kitimat Valley Airshed (TKVA) which stretchesfrom the small First Nations community of Hartley Bay in the south to Kalum Lake in thenorth. Aside from the locations listed above the airshed is sparsely populated. The airshed’s6complex topographical features along with the locations of its communities can be seen infigure 1.2.This valley is poised to undergo large-scale industrial expansion over the coming years; inaddition to the recently completed expansion of an existing aluminum smelter there is alsothe potential for the construction of:• a bitumen-condensate export-import terminal,• four natural gas liquefaction facilities, and• an oil refinery.All of these facilities emit air pollutants, many of which (especially NOX from liquefied naturalgas (LNG) processing and shipping), are O3 precursors. The current status of all proposedprojects is described in table 1.2.Table 1.2: Status of proposed projects in KitimatIndustry Project Status (August 2015)SmelterRio Tinto Alcan Kitimat permit approved, pending appeal andPrimary Works review by BC environmental appeal boardLiquefied Natural Gas*LNG Canada EA certificate granted June 2015Kitimat LNG EA certificate granted originally as LNG importfacility in October 2011, in permitting phaseDouglas Channel LNG EA Process (final decision October 2015)Other LNG (Triton LNG) EA ProcessOil RefineryKitimat Clean EA ProcessBitumen Import ExportEnbridge Gateway EA Certificate granted pending completionof 125 tasks*(British Columbia Ministry of Natural Gas Development, 2015b)1.2.1 Liquefied Natural GasAt -162 °C natural gas (mainly CH4) condenses into a liquid and occupies 1/600th of itsgaseous volume (Groupe International Des Importateurs de Gaz Naturel Liquifie, 2009). Onceliquefied, it can be loaded on to marine vessels and transported overseas for sale into export7Figure 1.2: The Terrace - Kitimat valley airshed. The green line outlines the Terrace - Kitimatvalley airshed and the thin brown lines are 250 m contours. ￿, ￿, ￿, ￿ and ￿ identifythe locations of Terrace, Thornhill, Kitimat, Kitamaat Village and Hartley Bay respectively.Parallel lines identify the locations of roads. Inset: The location of Kitimat within Canada.Grey lines show provincial and territorial borders. ￿ and￿ identify the locations of Kitimatand Vancouver respectively8markets. Natural gas is extracted primarily in northeastern BC in areas such as the HornRiver and Montney basins, and is proposed to be transported to the coast via pipelines (BritishColumbia Ministry of Natural Gas Development, 2015a).The provincial government claims that countries such as Japan and China have a thirst forLNG as replacement for nuclear and coal in power generation (Canadian Press, 2013); thisexplains in part why there are four proposed LNG projects in the TKV. Thus on the surfacethere seems to be some global benefit from utilizing this energy source. However, LNG pro-cessing has also raised alarm in and around communities where facilities are proposed becauseits process is emissions-intensive.Cooling natural gas to its condensation point requires tremendous energy; currently mostproponents are proposing to power their refrigeration processes by burning some of theirnatural gas in gas-fired turbines which drive compressors (Hoffman, 2015). Occurring at hightemperatures, natural gas combustion emits large quantities of NOX (for example, each LNGCANADA refrigeration circuit requires two turbines, each emitting NOX at a rate of 32.1 kg/hr(Reid et al., 2014)). If all proposed projects in the TKV are built, NOX emissions will reachapproximately 28 tonnes per day (refer to emissions chart - figure 2.11 in section 2.3.1, thisincludes the LNG facilities and their associated marine shipping emissions). To put this intocontext, correspondence with Environment Canada and Translink in June 2013 estimated thateach tonne of emitted NOX per day is equivalent to the NOX emissions from approximately22,000 vehicles in the Lower Fraser Valley (LFV) which includes the city of Vancouver and itssurrounding communities, each driving an average of 37 km per day (this includes a mix ofpassenger cars, passenger trucks and heavy duty trucks). Thus the total projected daily NOXemissions in the TKV are the equivalent of approximately 616,000 vehicles (Vingarzan, 2013).Assessing whether or not emissions of O3 precursor pollutants will actually lead to the produc-tion of O3 requires an air quality model capable of managing the complex chemistry presentedin section 1.1. These models are often referred to as chemical transport models or photochem-ical models, and are described briefly in the following section.91.3 Photochemical ModelsPhotochemical models are air quality models that simulate the production, loss and advectionof photochemical pollutant concentrations using a series of equations which characterize thechemical and physical processes in the atmosphere. Most of these models employ an Eulerian-based fixed coordinate system with (respect to the ground) consisting of a three-dimensionalseries of cells. At each time step in the model period, emissions advect from one cell to anotherand are subjected to the physical and chemical processes quantified in the model’s code (U.S.Environmental Protection Agency, 2015b). Final concentrations within each grid cell are afunction of six factors:1. emissions within the cell,2. deposition within the cell,3. chemical production in the cell,4. chemical loss inside the cell,5. advection of pollutants into the cell from adjacent cells, and6. advection of pollutants out of the cell into adjacent cells (Jacob, 1999).Chemical production and loss are highly nonlinear, as reaction rates for many compoundsdepend on environmental conditions as well as ambient concentrations themselves (Jacob,1999). A broad schematic of the photochemical grid modelling system is presented in figure1.3.1.3.1 Comprehensive Air Quality Model with ExtensionsThe Comprehensive Air Quality Model with Extensions (CAMx) is an Eulerian photochemicaland dispersion air quality model currently in its sixth edition. It is built and maintained byEnviron International Corp. (a global air quality consulting company) and is freely download-10Figure 1.3: A broad schematic of the Eulerian photochemical grid modelling system (Steynet al., 2012).able from the internet. Its FORTRAN code and input/output file formats are based on theUrban Airshed Model (UAM) convention. (ENVIRON, 2013). CAMx’s numerous data in-put requirements (many of which are also models) fall into five general categories: emissions,meteorology, photolysis, geographic and air quality. A schematic of the CAMx modellingframework including its pre and post-processors are presented in figure 1.4. Each category iscolour-coded. More information on the sub-models are presented in chapter 2. CAMx, withits Eulerian nature and chemical processes, is capable of assessing O3 and its precursors.11Figure 1.4: The overall modelling framework and input categories in CAMx (ENVIRON,2013).1.4 Research MotivationThe cumulative effects of airborne emissions from existing and potential sources in the TKVAwere the subject of a recent government study, however, despite high projected emissions ofphotochemical precursors, it was limited in scope to sulphur dioxide (SO2) and NO2 exposureas well as sulphur and nitrogen deposition. Primary particulate matter (PM) and secondarypollutants were omitted. Yet without an understanding of the total change in atmosphericpollutant concentrations, it is not possible to make final conclusions regarding either the risk(s)from any combination of new sources or the future status of the airshed.Research pertaining to photochemical modelling of O3 in sparsely populated areas with com-plex coastal terrain at high latitudes, particularly in the context of industrial development,is limited. A list presented in Simon et al. (2012) of 69 peer-reviewed articles concerning theapplication of photochemical models in North America between 2006 and 2012 suggests thatmost applications are regional or continental in nature, and only a handful are applied on a12local scale. Very few relate to industrial development and those that do are regional in scale.One similar study was conducted by Castell et al. (2010), who showed a range of potentialincreases to ambient O3 caused by a hypothetical natural gas power plant in a coastal Spanishvalley. Model output indicated O3 increases 20 - 40 km and 80 - 140 km away during oneepisode, and 20 - 60 km away during another. Spatially, O3 increases were greatest in NOX -sensitive areas (nearby forests and hillsides), i.e.: locations with emissions of high biogenicVOCs (bVOC) and low NOX emissions. Note that this study was conducted in a locationthat was neither sparsely populated nor situated at high latitudes, but did pertain to complexcoastal terrain and industrial development similar to those proposed for the TKVA.The most similar physical environment to the TKV where O3 has been extensively studied ona local scale is the LFV. The LFV is a temperate and triangular-shaped coastal valley boundedby the Coast mountains to the north, the Cascade mountains to the east and south and theGeorgia straight to the west (Armstrong, 1990). Steyn et al. (2013) and Ainslie et al. (2013)used a combination of models, observations and emission inventories to better understand therelationships between reductions in precursor emissions and episodic O3 concentrations. Modelsimulations indicated that O3 concentrations were highest in tributary valleys on the northernedge of the LFV in locations with no ambient O3 monitoring (Steyn et al., 2013). Furthermore,changes in valley-wide emissions over past decades have changed the O3 sensitivity of certaincommunities within the LFV. For example, the town of Chilliwack has generally changed frombeing VOC-limited to NOX -limited (Ainslie et al., 2013). Understanding this kind of changeis important as it has implications for airshed management.While the LFV is similar to the TKV in some respects they are also quite different. The TKVis much narrower and its terrain is more complex. As well, the population of the TKV is verysmall compared to the LFV (25,000 compared to approximately 2,500,000 (Statistics Canada,2015b)). Anthropogenic emissions in the TKV are confined to two small urban areas as wellas some mobile sources (roads and rail) while in contrast, anthropogenic emissions in the LFVoccupy a significant percent of the domain. Lastly, studies in the LFV have been broad in13nature reflecting large-scale changes to airshed emissions while questions in the TKV concernindustrial development from a handful of facilities and associated marine shipping emissions.Considering the lack of literature in this area there is currently an opportunity to apply aphotochemical model to study industrial development on a scale not often used and to aphysical environment that has received little attention. If successful, this approach can beapplied to other airsheds with complex coastal terrain subject to industrial expansion whereproposed projects emit photochemical precursor pollutants.1.5 Proposed Research QuestionsGiven the discussion in the previous sections, the following three research questions have beendeveloped:1. Can spring and summer photochemical O3 be replicated by a model in a manner thatis fit for the purpose of investigating worst-case concentrations that may result fromproposed industrial emissions in a constrained coastal airshed with complex terrain suchas the Terrace - Kitimat valley?2. Should all proposed industrial facilities in the Terrace - Kitimat valley airshed be con-structed,(a) Where are the locations of the O3 maxima, and what are the worst-case concentra-tions that could occur there?(b) What meteorological factors contribute to enhancing or reducing O3 concentrationsin the Terrace - Kitimat valley airshed?3. What is the current and future O3 sensitivity of the Terrace - Kitimat valley airshed toemissions of NOX and VOCs?Question 1 is broad in nature while Questions 2 and 3 are more specific. Successfully answeringQuestion 1 allows Questions 2 and 3 to be addressed. Question 1 involves the use of a control14case to demonstrate model fitness, i.e.: the modelling of known emissions and their levelof agreement with measured values. Question 2 involves the modelling of a test case withprojected emissions. The process of answering Question 2 includes not only an analysis offuture concentrations but also an investigation into the changes in pollutant concentrationsbetween control and test cases. This two step approach (whereby a model’s credibility inpredicting future states is based on its ability to describe a current state) is not uncommon(Song et al., 2010). Question 3 will increase our understanding of the chemical sensitivity ofthe valley, and may inform O3 reduction strategies in future airshed management planninginitiatives.15Chapter 2MethodsRunning CAMx in order to answer Questions 1, 2 and 3 requires up-front groundwork, cate-gorized into three groups:1. Period selection.Period selection is the process of narrowing the model simulation to only those timesconsidered necessary to address the research questions being asked, as computationalresources are expensive. Conditions leading to worst-case O3 must be identified basedon some predictors and model runs must be constructed around those times.2. Developing and evaluating input files.Developing and evaluating input files is the process of creating and evaluating all of theinputs used by CAMx. As displayed in figure 1.4, the CAMx modelling system relieson a number of inputs which must be carefully developed in order to obtain reasonableresults. Developing inputs requires at least as much time as running CAMx itself asmany of these inputs are also models. Output from these models can only be used asinput to CAMx after some form of fitness evaluation; in this case some outputs wereevaluated quantitatively while others were evaluated qualitatively.163. Model settings.Model settings are the various switches and parametrization schemes in CAMx that caninfluence the model’s outcomes. For example, the CB05 chemical mechanism in CAMxcontains 51 chemical species and 156 reactions while the CB6 mechanism contains 77 gasphase species and 218 reactions (Yarwood et al., 2010). Professional experience oftenguides the selection of one scheme over another as well as the nature of the emissionsinventory and research questions (UNC, 2013).Completing this up-front groundwork comprises the modelling approach. This chapter de-scribes the model approach for each group listed above.2.1 Period SelectionTwo periods were selected for this research, one in the spring and the other in the summer.A spring period was chosen because O3 is naturally elevated during this time (Monks, 2000)and the result of adding additional NOX into the TKVA warranted some investigation. Asummer period was selected as summer is a more traditional time of year for O3 episodes dueto increased production of the OH radical through reactions 1.3 and 1.4 (Jacob, 1999).Figure 2.1 illustrates the annual O3 cycle in Smithers, BC, the closest location to the TKVwith ambient O3 monitoring. The mean maximum daily O3 mixing ratios for the years 2007through 2013 along with daily maximum values form three select years, 2008, 2011 and 2012are presented. The springtime O3 peak is evident in this figure. Note that there is no summerpeak in Smithers and it is probable that communities in the TKVA share this annual O3 cycle.The year 2010 was selected for this research due to the availability of emissions data from theRio Tinto Alcan (RTA) aluminum smelter and the 2010 National Pollutant Release Inventory(NPRI) national emissions inventory. Finding a range of dates for each period was desiredwhere either:(a) O3 was already elevated, or17Figure 2.1: Maximum daily O3 mixing ratios at Smithers. Black line is the mean seven yearperiod (2007-2013), red dotted line is 2008, blue dotted line is 2011 and orange dotted line is2012.(b) conditions conducive for photochemical O3 production were present.As no lengthy O3 dataset exists for the TKV (permanent O3 monitoring only began in Terracein the spring of 2015), exact dates for the model period were informed by a random forests (RF)regression model which incorporated measured ambient meteorological and air quality variablesas O3 predictors. The model was trained with data from Smithers and was then applied toKitimat. Note that while the climates of the TKV and Smithers are different (coastal andcontinental respectively), research suggests that the general O3 climatology is similar across thenorthern hemisphere, particularly the presence of a springtime peak (Monks, 2000). A reviewof the limited 2015 springtime O3 data in Terrace confirms the presence of the springtime peakwhile a review of summertime O3 data in Kitimat (from a temporary monitoring program)suggests that afternoon peak mixing ratios are similar to Smithers, while overnight titrationis not as prevalent (not shown).A complete description of this process is presented in appendix A, where it is shown that themodel estimated elevated O3 in 2010 during two periods in the spring (April 14th - 18th andMay 2nd - 14th), and two periods in the summer (July 29th - Aug 6th and Aug 12th - 18th).These days coincided with sunny weather and low relative humidity (RH). The final model18period for each season was selected to start at the beginning of the first period of the seasonand continue to the end of the second. The dates between were included for conveniencehowever were not expected to be conducive to O3 production. Table 2.1 shows the final datesfor both the spring and summer periods.Table 2.1: Dates for the spring and summer model periods (2010).Season Start Date End Date Number of DaysSpring April 14 May 14 31Summer July 29 August 18 212.2 Developing Meteorological InputsAs can be seen in figure 1.4, successfully running all other components of the photochemicalmodelling system requires completed meteorological files as inputs in some manner (with theexception of the photolysis component). Meteorological inputs are perhaps the most importantinputs to any dispersion or chemical model as without reliable meteorological inputs they willfail to produce realistic results. This section concerns information related to the meteorologicalinputs used in this research.2.2.1 Weather Research and Forecast ModelExcluding data collection (from meteorological monitoring instruments, satellites, etc.), thefirst step in generating meteorological inputs for any photochemical model is usually to runa prognostic weather model. For this research, meteorological inputs were generated usingthe Weather Research and Forecast Model (WRF) model. WRF is a dynamical mesoscaleprognostic numerical weather prediction model developed in an ongoing collaboration in-volving the National Centre for Atmospheric Research (NCAR) and the National Oceanicand Atmospheric Administration (NOAA)’s National Centre for Environmental Prediction(NCEP), with several other US government agencies and universities and its user community(National Oceanic & Atmospheric Administration, 2014). According to the WRF website19(http://www.wrf-model.org/index.php), WRF has over 25,000 users globally in more than130 countries.For the purpose of this research WRF (ARW core, version 3.6) output was obtained from DavidSuita, a PhD student working at the Geophysical Disaster Computational Fluid DynamicsCentre at UBC. Each run consisted of a 34-hour run, initialized at 00:00 UTC (16:00 PacificStandard Time (PST)). Output from the first eight hours was discarded as model spin-up,and the final two hours were also removed, leaving 31 (spring) and 21 (summer) 24-hourperiods. (CAMx modelled each day separately and used final conditions from one day asinput conditions for the following day.) WRF was initialized with 32 km North AmericanMesoscale Forecast System (NAM) gridded output. WRF output was developed for a parentand series of nested grids using a Lambert Conformal projection with horizontal resolutionsof 36, 12, 4 and 1.333 km. Figure 2.2 shows the locations of the four grids and their positionsrelative to North America.Each domain has a number of east-west and north-south cells (x and y-cells respectively), aswell as vertical layers (z-cells). Table 2.2 lists the properties of each model domain.Table 2.2: Properties of the WRF model domainsDomain Name Horiz. resolution x-cells y-cells z-cells Grid centre lat Grid centre lonParent d01 36 km 99 99 40 54.230 -128.6201st nested d02 12 km 99 99 40 53.883 -129.1572nd nested d03 04 km 120 105 40 54.012 -128.5093rd nested d04 1.333 km 120 159 40 54.117 -128.706The smallest domain, d04, stretches 206.7 km in the north-south direction and 156 km in theeast-west direction. This orientation reflects the orientation of the TKVA. Figure 2.3 shows anoutline of d04 along with the locations of TKVA communities Terrace, Kitimat, and HartleyBay.From a north-south perspective, Kitimat is located in the middle of the domain while Terraceis located in the centre of the northern half. This was done intentionally to locate the majorityof the emissions in the centre of the domain where they could advect north or south dependingon wind direction. Hartley Bay, the small First Nations community at the southern edge of20Figure 2.2: All WRF domains. Green and red lines identify the 36 km and 12 km domainsrespectively, while the blue and black rectangles identify the 4 km and 1.333 km domainsrespectively.the TKVA, is located near the southern end of the grid.Model SettingsSome notable settings for the WRF model physics are provided in table 2.3. These wereselected based on David Suita’s judgement and experience.21Figure 2.3: Outline and individual cell structure of the WRF d04. The maroon line is thedomain outline and the thin black lines are the 1.333 km cells. ￿, ￿ and ￿ identify thelocations of Terrace, Kitimat and Hartley Bay respectively. Parallel lines identify the locationsof roads.Table 2.3: Physics options used in WRFPhysics Setting WRF name Option Scheme nameLongwave radiation ra_lw_physics 1 RRTMShortwave radiation ra_sw_physics 1 MM5 shortwaveSurface layer sf_sfclay_physics 1 Monin-Obukhov similarity theoryLand/water surface sf_surface_physics 2 Noah land surface modelPBL physics bl_pbl_physics 1 YSU PBL schemeCumulus parameterization cu_physics 1 New Kain-FritschMicrophysics mp_physics 4 WSM 5-classTurbulence / Diffusion diff_opt, km_opt 1, 4 2nd order diffusion on model levels222.2.2 Weather Research and Forecast Model EvaluationWRF output was evaluated both quantitatively and qualitatively. The quantitative evaluationinvolved a comparison of modelled and observed parameters at three locations in the TKV,while the qualitative evaluation judged the spatiotemporal evolution of model output over thecourse of a single day. The complete evaluation can be found in Appendix C, where it isdemonstrated that this WRF output is fit for use in CAMx to address the research questionsfrom section 1.5. A summary of the evaluation methodology is presented in the followingsection.Quantitative EvaluationFor the quantitative evaluation, modelled vs. observed variables were compared at threemeteorological monitoring stations in the TKV. Only output from d04 was used. Propertiesof these stations are listed in table 2.4. Note that for all stations the actual elevation wasgreater than the elevation assigned to the station’s grid cell in WRF.Table 2.4: Properties of the meteorological monitoring stations used for the WRF evaluation.Location Lat Lon Elevation Elevation (WRF) x-cell y-cellm-asl m-asl WRF WRFTerrace Access Centre 54.518 -128.598 68 62 66 111CYXT Terrace Airport 54.471 -128.573 210 170 67 107Kitimat Whitesail 54.069 -128.639 92 83 63 75The locations of these stations are shown in figure 2.4. While the Terrace B.C. Access Centreand CYXT airport are not far each other (only 5 km), the difference in elevation betweenstations as well as the position of Terrace relative to the Skeena river valley led to verydifferent results in the quantitative evaluation.The parameters evaluated were temperature, RH, pressure, wind direction and wind speed.Many of these were proposed by Dennis et al. (2010) as being important from a meteorologicalevaluation perspective and have the advantage of having been measured at multiple locations23Figure 2.4: Evaluation locations for WRF d04 output. ￿, (, and ￿ identify the locationsof Terrace BC Access Centre, CYXT and Kitimat Whitesail respectively. The maroon lineoutlines the WRF d04 and the green line outlines the Terrace - Kitimat valley airshed. Parallellines identify the locations of roads.across the valley. Parameters measured at each station and used in the evaluation are listedin table 2.5.24Table 2.5: Parameters measured at meteorological monitoring stations in the Terrace - Kitimatvalley airshedLocation Temp RH Pressure Wind Speed Wind Direction°C % Pa ms-1 degTerrace Access Centre X X X XCYXT Terrace Airport X X X X XKitimat Whitesail X X XAs demonstrated in appendix C, the quantitative evaluation of WRF output shows that themodel’s performance was fit for use in CAMx. Results were strongest for pressure and tem-perature and weakest for RH. Results for wind direction and wind speed were mixed. Despiteoverall strength of the model temperature output, one important limitation worth noting wasthe model’s inability to recreate peak afternoon temperatures. On days where the measuredtemperature was grater than 30 °C, WRF was unable to recreate the highs of the diurnalcycles, typically falling short by 5 °C. As discussed later in section 4.2.3, this creates someuncertainty with respect to VOC emissions and resulting O3 mixing ratios, as emissions ofbVOC increase exponentially with increasing temperature, particularly isoprene (Guentheret al., 1993).Qualitative EvaluationThe qualitative evaluation consisted of a spatiotemporal evaluation of wind streamlines andtemperature evolution over a 24-hour period. This evaluation was based on personal experiencewith the TKV and flow in complex terrain. The day of July 31st was selected as this was a hotsunny day with a measured diurnal offshore-onshore wind. The qualitative evaluation showedthat both wind streamlines and temperature evolved intuitively over the course of space andtime for the chosen day in the evaluation. Streamlines were concentrated in overnight hoursand there was evidence of katabatic flow on one side of the TKV. The offshore breeze shifted tobecome an onshore breeze around noon, about the time one would expect for a day in July inKitimat, and the model was able to resolve flow around very fine-scale topographical featureslike tertiary valleys. The spatial and temporal evolution of temperature was as expected for a25summer day in complex coastal terrain, where the land 2 m temperature was originally coolerthan the 2 m temperature of the nearby ocean though this reversed around noon. Temperaturesat higher elevations also increased to some extent, albeit with a smaller amplitude than thevalley bottom temperature.2.2.3 Intermediate Programs: MCIP and WRFCAMxThe second step in developing meteorological inputs for a photochemical model involves trans-lating the relevant outputs from the prognostic weather model into a format used by thephotochemical model. This is done by taking the meteorological model output fields in theirnative formats and: performing horizontal and vertical coordinate transformations, diagnosingadditional atmospheric fields, and defining gridding parameters (Otte and Pleim, 2010). TheMeteorology-Chemistry Interface Processor (MCIP) (as described in Otte and Pleim (2010))and WRFCAMx (as described in ENVIRON (2013)) fulfill this role. As some of the inputs toCAMx were built for use in Community Multiscale Air Quality Modelling System (CMAQ),CAMx requires both MCIP and WRFCAMx to be used. Table 2.6 sorts the components ofthe CAMx modelling system according to the translator program used (refer to 1.4 to viewthese in the overall context of CAMx).Table 2.6: CAMx input models sorted by translator programWRFCAMx MCIPCAMx SMOKE(anthropogenic emissions model)MOZART MEGAN(boundary and initial (biogenic emissions model)conditions model)MERGE_LULAI(merges land use andleaf area index files)SEASALT(marine emissions (mainlyaerosol phosphate) notused in this research)In addition to the functions described above, MCIP and WRFCAMx also act to shrink the26WRF domain into the CAMx domain in both the horizontal and vertical directions. TypicallyCAMx domains are between 4 and 12 cells narrower (i.e.: between 2 to 6 cells per side) thanWRF domains and contain fewer vertical levels. Table 2.7 lists the properties of the CAMxmodel domains defined in MCIP and WRFCAMx. An outline of the 1.333 km CAMx domainrelative to the WRF domain can be seen in figure 2.5, while figure 2.6 displays the mappingof vertical layers from WRF to CAMx.Table 2.7: Properties of the CAMx model domainsDomain Name Resolution x-cells y-cells Vertical Layers Grid centre lat Grid centre lonParent grd01 12 km 93 93 23 53.883 -129.1571st nested grd02 04 km 116 101 23 54.012 -128.5092nd nested grd03 1.333 km 110 155 23 54.117 -128.706Note that the 36 km domain was not included in CAMx. It was needed in order to developthe WRF input files (i.e.: to ensure large scale atmospheric processes trickled down to thefiner grids) but was not required for the local-scale chemical modelling.Also note that the 1.333 km CAMx domain is ringed by one extra set of cells (not shown),referred to as buffer cells. These cells hold internal boundary conditions between a finer gridand the coarser grid surrounding it. (The 4 km CAMx domain also contains a ring of buffercells.)As can be seen from figure 2.6, the transition from 40 to 23 vertical layers occurs gradually.The lower 10 layers of the WRFCAMx domain is structured identically to the lower 10 WRFlayers. From layer 11 to 17, WRFCAMx layers occupy two WRF layers and above thisWRFCAMx layers occupy three WRF layers (with the exception of the highest WRFCAMxlayer which only occupies two WRF layers). This layer structuring was determined based oncorrespondence with ENVIRON staff (Wilson, 2014).27Figure 2.5: Outlines of the WRF d04 and CAMx grd03. The maroon line outlines WRFd04 and the steel blue line outlines the CAMx grd03. The green line outlines the Terrace -Kitimat valley airshed. ￿, ￿ and ￿ identify the locations of Terrace, Kitimat and HartleyBay respectively. Parallel lines identify the locations of roads.28Figure 2.6: The vertical levels in WRF mapped to the vertical levels in CAMx. The ordinateof a) shows the ETA values of the WRF while the ordinate of b) shows the ETA values ofthe CAMx. The abscissa in both figures shows the number of vertical levels in the respectivemodel. From this figure it is possible to see, for example, that the 11th layer in the CAMxcorresponds to the 11th and 12th layers in the WRF. Respective layer ETA midpoints aredenoted with dots.2.3 Developing Emissions InputsOnce meteorological files are developed it is possible to move onto the next important set ofinputs, emissions. As these research questions revolve around the effect of proposed emissionsit is necessary to ensure they are adequately characterized. Two emissions models were usedto calculate and process emissions for the TKVA airshed:• The Sparse Matrix Operator Kernel Emissions Modelling System (SMOKE) and• the Model of Emissions of Gases and Aerosols from Nature (MEGAN).More information on both of these models is given in the following sections.2.3.1 Anthropogenic EmissionsAnthropogenic emissions were processed in SMOKE, an emissions processor. SMOKE wascreated in order to prepare emission inventories for use in air quality and chemical models,29which involves transforming an inventory through temporal and spatial allocation as wellas chemical speciation (UNC, 2013). SMOKE was created at the Microelectronics Centreof North Carolina’s Environmental Monitoring Centre and is now distributed through theCommunity Modelling and Analysis System (CAMS) operation (https://www.cmascenter.org/) at Chapel Hill in North Carolina. SMOKE processes all criteria air contaminants(CAC) as well as a variety of toxics (benzene, formaldehyde, etc.). It supports a variety ofanthropogenic emission source types such as area, point, mobile as well as biogenic sources.The following subsections concern information regarding the temporal, spatial and chemicalspeciation of both the Control (used in addressing Questions 1, 2 and 3), and Test (used inaddressing Questions 2 and 3) cases.Control CaseThe Control case (also referred to as ’Base’) served two functions and was used for:1. Model evaluation, and2. Calculating percent change between Control and Test CasesIn the TKVA, where the Control case was needed for the second of the above functions,emissions from three source-types were modelled:1. industrial sources2. mobile sources, and3. biogenic sources.The industrial sources consisted of a single permitted facility, the RTA aluminum smelter,located to the southwest of Kitimat. Within this one facility were many point sources as wellas some line sources (treated in SMOKE as area sources based on correspondence with otherprofessionals (Henolson, 2014)). The location of anthropogenic sources can be seen in figure2.7. In this image roads and rail lines should be considered as emission sources.30Figure 2.7: Emission locations in the Terrace - Kitimat valley airshed for the Control case.identifies locations of aluminum smelter sources while parallel lines and the solid black lineindicate locations of mobile sources (roads and rail sources respectively). ￿,￿ and￿ identifythe locations of Terrace, Kitimat and Hartley Bay respectively. The green line outlines theTerrace - Kitimat valley airshed. Inset: A closer look at emission locations in Kitimat.31Emissions information for RTA were obtained largely from the BC Ministry of Environment(MOE) (ESSA-Technologies et al., 2014). CO emissions (for some sources) were obtainedfrom the NPRI, while VOC emissions were not available from either BC MOE or the NPRI.Omission of the smelter’s VOC emissions adds uncertainty to the Control case model scenario.Mobile sources consisting of road (highway only) and rail emissions were derived from the 2010SMOKE-ready emissions inventory produced by the NPRI, available at http://www.epa.gov/ttnchie1/net/canada.html. Mobile emission totals were available for the Skeena region ofBC and were allotted proportionately to cells with roads or rail lines based on the fraction ofroad or rail in the cell with the total length of roads and rail lines in the region. Daily CACemissions from the above sources are presented in figure 2.8. (Note: NOX emissions are oftenapportioned a NO-NO2 split of 90%-10% (Swedish Environmental Protection Agency, 2009).)Figure 2.8: Daily criteria air contaminant emissions in the Terrace - Kitimat valley airshedfor the control case.As there was insufficient O3 monitoring in either Kitimat or Terrace, model evaluation occurredin Smithers, the closest location with ambient O3 data. Smithers is located in the BulkleyValley (BV), adjacent to the TKV and separated by the Coast Mountain range. The BV con-sists of five incorporated communities, the largest of which is Smithers, pop. 5400 (StatisticsCanada, 2014b). Similar to the TKV, emissions from three source-types we modelled as part32of the model evaluation:1. industrial sources2. mobile sources, and3. biogenic sourcesThe industrial sources consisted of two facilities within Smithers and two to the southeastin Houston (65 km away). Of the facilities in Smithers one was a sawmill and the other aparticleboard plant. Sources in Houston included a sawmill and a pellet plant. The locationof BV anthropogenic sources can be seen in figure 2.9.Point source emissions data were obtained from the 2010 SMOKE-ready emissions inventoryproduced by the NPRI. Possibly due to a lack of reporting requirements, emissions from thepermitted sources had disparities with those obtained from the BC MOE permit fee database(Pierce, 2015). Furthermore, these industrial point sources contained no stack information(height, exit velocity, exit temperature, etc.) and so were eventually treated by CAMx asarea sources. Despite these differences, the NPRI emissions were used to maintain a level ofconsistency in the approach (i.e.: without stack testing for NOX it is hard to know whichnumbers to believe). The BV mobile source emissions inventory was developed in the samemanner as for the TKVA. Daily CAC emissions form the above sources are presented in figure2.10. Note that the point source emissions from Houston are omitted in this inventory.Test CaseThe purpose of the Test case was to investigate: worst-case O3 concentrations, the locationsof those maxima, the meteorological factors that contributed to enhancing or reducing O3,and the future O3 sensitivity of the TKVA to emissions of NOX and VOCs. The Test caseconsisted of the same source types as the Control case (industrial, mobile and biogenic) thoughadditional sources were included:• an expanded aluminum smelter (RTA),33Figure 2.9: Outlines of the Bulkley and Terrace - Kitimat valley airsheds and emission sourcelocations in the Bulkley valley for the Control case. The thick light blue line outlines theBulkley valley arished while the thick red line outlines the Terrace - Kitimat valley airshed.identifies locations of emission sources in Houston and identifies the location of an emissionsources in Smithers (there are two sources but their locations overlap in the larger map). ￿,￿,￿,￿ and￿ indentify the locations of Terrace, Kitimat, Hazelton, Smithers and Houstonrespectively. Parallel lines identify the locations of roads while the solid black line identifiesthe locations of railways. Inset: emission source locations in Smithers for the control case.identifies the location of the sawmill while identifies the location of the particleboard plant.Parallel lines identifies the locations of roads while the solid black line identifies the locationsof railways. ￿ identifies the location of the St. Josephs air quality monitoring station.34Figure 2.10: Daily criteria air contaminant emissions in the Bulkley valley airshed for thecontrol case.• four LNG facilities including:1. LNG Canada,2. Kitimat LNG,3. Douglas Channel LNG, and4. Triton LNG,• an oil refinery, and• marine emissions associated with shipping from all sources.All industrial sources were input as point sources (SMOKE transforms some point sourcesto area sources in later processing). Information for the additional emssions was obtainedfor a variety of sources. For LNG Canada, emissions information was obtained from theirEnvironmental Assessment application (Reid et al., 2014). For all other point sources emis-sions information was obtained from the BC MOE (ESSA-Technologies et al., 2014), with35the exception of CO and VOC emissions. These data had to be derived based on scaling ofequivalent sources from the LNG Canada application, as that inventory was more complete.In some cases information from the NPRI was used to supplement missing information. If itwas not possible to derive VOC or CO emissions they were set to 0. To simulate ship move-ment marine emissions were considered as a series of 27 stationary points within the Douglaschannel (this same approach was used by ESSA-Technologies et al. (2014)). The locations ofthe new sources are not confined to the community of Kitimat and are spread up and downthe TKV. The northernmost source is the oil refinery while the marine emissions are spreadacross the Douglas channel, reaching south past Hartley Bay. These locations are presentedin figure 2.11.Daily CAC emissions for the Test case are presented in figure 2.12. Emissions of all CACsincrease relative to the Control case, some by over 700%. Of all the potential new sources tothe TKVA, the largest emitters of the photochemical precursor NOX are the marine sourcesthat transport LNG. As mentioned, marine sources are spread over the length of the Douglaschannel and each individual source is 1/27th of the total emissions. Also recall that theseemissions represent the combined emissions from all projects, as the marine emission inventorywas developed in this manner by ESSA-Technologies et al. (2014) for the BC MOE. The projectwith the highest NOX emissions is the LNG Canada project.2.3.2 Biogenic EmissionsBiogenic Emissions in the TKVABiogenic emissions, particularly bVOC, play a critical role in the photochemistry of the TKVAand indeed North America. In the USA annual bVOC emissions exceeded anthropogenic VOCemissions (Jacob, 1999). In the TKVA biogenic emissions are by far the largest source ofVOCs (by two orders of magnitude). Isoprene, the largest bVOC emitted in the TKVA, isvery reactive, with a lifetime of 0.2 days (Fall, 1999) and has been found to produce O3 at0.182 ppb/hr (rate coefficient against OH, kOH x 1012 = 101) in England (Derwent, 1999).36Figure 2.11: Emission locations in the Terrace - Kitimat valley airshed for the Test case.identifies locations of aluminum smelter sources, identifies locations of liquefied natural gassources, identifies locations of refinery sources and identifies locations of marine sources.Parallel lines and the solid black line indicate locations of mobile sources (roads and railsources respectively). ￿,￿ and ￿ identify the locations of Terrace, Kitimat and HartleyBay respectively. The green line outlines the Terrace - Kitimat valley airshed. Inset: A closerlook at emission locations in Kitimat.37Figure 2.12: Daily criteria air contaminant emissions in the Terrace - Kitimat valley airshedfor the Test case.In order to appreciate the importance of biogenic emissions it is necessary to understand thefactors that contribute to natural emissions, namely, the forest and land cover types in theTKVA. Forest type and landcover in the TKVA are themselves functions of the biogeoclimaticzones in which the valley is located.The TKV bottom is located within the Coastal Western Hemlock Biogeoclimatic EcosystemClassification zone of BC (Meidinger and Pojar, 1991), (Banner et al., 1993). As elevationincreases there is a transition to the Mountain Hemlock Biogeoclamatic zone while at highelevations (above 1000 m) the Alpine Tundra zone is prevalent. These zones are depicted infigure 2.13 (a). The vegetation in this area ranges from coastal and valley-bottom wetlands andbog-forests, to floodplain shrubs and deciduous forests, to coniferous forests with understoriesthat typically include well-developed shrub and moss layers in stands with open canopies, andfinally subalpine scrublands. Forest communities vary with elevation and moisture availability,however overall forest productivity is high in the region. Low-lying coastal areas tend toinclude western red and yellow-cedar, western hemlock and Sitka spruce (ESWG, 1995). Forestharvesting has converted much of the valley bottom (approximately 97% (Douglas, 2012)) to20-50 year old second growth western hemlock regeneration that typically creates a very dense38canopy and an impoverished understory. Floodplains include dense shrub cover of red-osierdogwood, red elderberry, salmonberry, horsetails and ferns (Demarchi, 2010). Red alder andblack cottonwood are common in disturbed alluvial landforms. Western hemlock occurs upto approximately 460 m elevation with understory species such as blueberry, false azalea anddevil’s club (ESSA-Technologies et al., 2013). At higher elevations a subalpine zone of yellowcedar, mountain hemlock and amabilis fir with an understory of blueberry and feather mossesare present. Wet meadows can also be found on perched benches that intercept groundwater.These species transition to stunted trees known as "krummholz", and beyond 1000 m elevation,to mountain heather hearthland and treeless alpine tundra (ESWG, 1995). Rock, snow andice, with pockets of alpine vegetation characterize the highest elevations (ESSA-Technologieset al., 2013). The various land cover types in the TKV are shown in figure 2.13 (b).39Figure 2.13: Classification of the Terrace - Kitimat valley for: a) biogeoclimatic zones and b) landcover types. ￿ and ￿ identifythe locations of Terrace and Kitimat respectively.40Table 2.8: Biogenic volatile organic compound species output by the Model of Gases andAerosols from NatureAbbreviation Full Group NameALD2 AcetaldehydeALDX Higher aldehydesETH EtheneEOTH EthanolFORM FormaldehydeIOLE Internal olefin carbon bondISOP IsopreneMEOH MethanolNR Nonreactive VOCsOLE Terminal olefin carbon bondPAR Parrafin carbon bondTERP TerpeneTOL Toluene and other monoalkyl aromaticsXYL Xylene and other polyalkyl aromaticsMEGANBiogenic emissions were developed and processed using MEGAN, a global biogenic emissionsmodel with a base resolution of approximately 1 km (WSU, 2012). The MEGAN user’s manualoffers an excellent description of its operation:MEGAN is a semi-mechanistic model that accounts for the major known processescontrolling biogenic emissions... Emissions of 150 chemical species are includedin MEGANv2.1 and the model can output individual compounds or categoriesassociated with various atmospheric chemistry schemes. The 150 compounds arelumped into 20 categories based on how emissions vary in response to changes inenvironmental conditions. Emission variations are first estimated for the 20 cate-gories and then speciated into the 150 compounds or output in chemical categoriesassociated with common atmospheric chemistry schemes (e.g., CB4, CB05, CB6,SAPRC99, MOZART, SOAX). Driving variables include land cover, weather, andatmospheric chemical composition. -MEGAN User’s Guide (WSU, 2012)Daily airshed-wide biogenic emissions created by MEGAN are shown in figures 2.14 (a) and2.14 (b) for the spring and summer model periods respectively. Table 2.8 lists the abbreviationsand full compound grouping for each bVOC output by MEGAN.41Figure 2.14: Daily biogenic volatile organic compound emissions in the Terrace - Kitimatvalley airshed for: a) spring and b) summer model periods.As can be seen from figures 2.14 (a) and 2.14 (b), summer bVOC emissions are approximatelyone order of magnitude greater than spring emissions. Airshed-wide bVOC emissions increasewith increasing temperature and exceed 300 tonnes per day when the temperature surpasses 29°C. Emissions are dominated by isoprene and other multiples of the C5H8 chain (e.g.: terpenes).Isoprene emissions come largely from deciduous forest stands (aspen trees are particularly highemitters (Sharkey et al., 2008)) and shrubs, and isoprene is thought to be synthesized andemitted as protection against heat flecks, very rapid changes in leaf temperature caused bysunlight (Sharkey et al., 2008). Figure 2.15 shows that biogenic isoprene emissions occur in42high quantities on the west side of the TKV. Cross-referencing this figure with figure 2.13 (b),it is possible to confirm that isoprene emissions are highest in areas with broadleaf, wetlandand scrubland land cover. Mature coniferous forest stands emit little isoprene in the TKVA,though there are high emissions on the southern edge of the domain which is designated asconiferous forest.Biogenic Emissions in the BVSimilar to the TKV, biogenic emissions play a critical role in the photochemistry of the BV.The BV is comprised of four distinct biogeoclimatic zones. The Alpine Tundra zone is above1500 m elevation while between 1000 m and 1500 m is the Engelmann Spruce - Subalpine Firzone (Meidinger and Pojar, 1991). Below 1000 m elevation two zones split the valley withthe northwestern quarter belonging to the Interior-Cedar Hemlock zone and the remainderbelonging to the Sub-Boreal Spruce zone (Banner et al., 1993). The Sub-Boreal Spruce zoneis characterized by cold winters and warm summers, deep snow cover and dense forests (Mei-dinger and Pojar, 1991). Coniferous tree species include hybrid white spruce, subalpine firand lodgepole pines, while deciduous species include trembling aspen, black cottonwood andpaper birch. Low elevation wetlands are common in this zone and consist of sedges, scrubbirch and willows, with white and black spruce on treed sites. Marshes occur around lakesand streams, usually with horsetails, sedges, cattail and bulrush (British Columbia Ministryof Forests, 1998c). The ecology of the Interior-Cedar Hemlock zone is notably different andis well defined in the BV by the presence of western hemlock. Deciduous tree species includetrembling aspen, and paper birch in dry areas. In the rough, steep terrain above these zones,forests of Englemann spruce and subalpine fir comprise the Engelmann Spruce - SubalpineFir zone. In the BV this begins approximately above 1000 m elevation. The zone containsfew deciduous trees of any kind. Shrubs such as huckleberry and grouse berry are present inopen areas and meadows (British Columbia Ministry of Forests, 1998b). Finally, the AlpineTundra zone is at the highest elevations. This is a treeless environment with steep and rockyterrain and snow-capped peaks. Vegetation is scarce and occurs largely in the mid and lower43Figure 2.15: Isoprene emissions in the Terrace - Kitimat valley on August 15th, 2010 at 16:00PST. The locations of Terrace and Kitimat are identified with ￿, and ￿ respectively.elevations of the zone. Plants of the Alpine Tundra zone are small, close to the ground andsnow-covered for up to 10 months of the year (British Columbia Ministry of Forests, 1998a).The delineation of these zones, along with the biogeoclimatic zones of the TKV are shown44in figure 2.16 (a). Lower elevations of the BV support extensive agricultural development,primarily cattle grazing and hay crops, resulting in land cover classification that differs con-siderably from the TKV where commercial agriculture is a minor land use component. Thiscan be seen by comparing 2.16 (b) with 2.13 (b). Once prevalent throughout, forest harvest-ing in the BV is now concentrated in side-valleys as agricultural lands have replaced a largeproportion of the historically forested lower elevations.45Figure 2.16: Outlines of the Bulkley and Terrace - Kitimat valley airsheds and classification of the region for: a) biogeoclimaticzones and b) landcover types (Bulkley Valley only). ￿, ￿, ￿, ￿ and ￿ indentify the locations of Terrace, Kitimat, Hazelton,Smithers and Houston respectively. Parallel lines identify the locations of roads.46Biogenic emissions developed and processed using MEGAN for the BV are presented in figures2.17 (a) and 2.17 (b) for the spring and summer periods respectively. In both figures theordinate has been scaled similarly to figures 2.14 (a) and 2.14 (b) for easy comparison withthe TKV. Based on land cover, higher valley bottom elevation (i.e.: cooler temperatures) anda smaller geographic area, bVOC emissions in the BV are approximately 50% of those in theTKV.Figure 2.17: Daily biogenic volatile organic compound emissions in the Bulkley valley airshedfor: a) spring and b) summer model periods.To emphasize the difference in biogenic emissions between the TKV and the BV, figure 2.18shows isoprene emissions in both valleys for the same hour as figure 2.15. From this image47it is apparent that the land cover in the BV does not produce the same level of isopreneas the TKV. Indeed, Sharkey et al. (2008) explain that crop plants which are selected forrapid growth require open stomata, and that high stomatal conductance allows high rates oflatent heat loss, buffering against heat flecks. Thus, crop plants should not, and generallydo not emit isoprene. Also, a review of Smithers and Kitimat summertime temperatures(not shown) suggests that Kitimat experiences higher temperatures than Smithers on dayswhere the maximum temperature is greater than 30 °C. Given that isoprene emissions increaseexponentially with increasing temperature (for more see section 4.2.3), on hot days one wouldexpect higher isoprene emissions from isoprene-emitting vegetation in the TKV than the BV.Note that in this image, the emissions scale has been adjusted to accommodate the increasein grid size for the 4 km domain (it takes 9 - 1.333 km cells to make 4 km cells so the emissionrate has been scaled to 0 - 900 mol/hr from 0 - 100 mol/hr).2.3.3 Other Emissions ProcessingIn sequential order, biogenic emissions were processed before anthropogenic emissions. Thisis because the final step of the SMOKE modelling system is the merging of anthropogenicand biogenic emissions. All non-point source emissions are output in a CAMx-usable format,and the final step in the emissions processing is to run the point source output files through aseparate processor, the PiGSET program. The modeller uses this program to select and set thesources which will be treated with the Plume-in-Grid (PiG) submodel in CAMx. PiGSET alsoconverts ASCII point source files commonly generated by emission models such as SMOKE toCAMx-ready binary format (ENVIRON, 2013). Once PiGSET has completed processing thepoint sources, all emissions output is ready to be used as input to the photochemical model.(Note that a discussion on plume-in-grid modelling occurs in section 2.6.4.)48Figure 2.18: Isoprene emissions in the Bulkley and Terrace - Kitimat valleys on August 15th,2010 at 16:00 PST, and outlines of the Bulkley and Terrace - Kitimat valley airsheds. Thethick light blue line outlines the Bulkley valley arished while the thick red line outlines theTerrace - Kitimat valley airshed. ￿, ￿, ￿, ￿ and ￿ identify the locations of Terrace,Kitimat, Hazelton, Smithers and Houston respectively.492.4 Developing Boundary and Initial Conditions InputsBoundary and Initial conditions (ICs) were generated from the Model for Ozone and Re-lated Chemical Tracers (MOZART), a global chemical transport model independent of otherchemical models. Detailed information about MOZART is published in Emmons et al. (2010).MOZART output includes 85 gas-phase species and 12 bulk aerosol compounds. 6-hour outputcan be downloaded for free from http://www.acom.ucar.edu/wrf-chem/mozart.shtml andcustomized for any domain. MOZART includes 39 photolysis and 157 gas phase reactions in itschemical process. The freely available output uses anthropogenic emissions based on the ARC-TAS global emissions inventory (http://bio.cgrer.uiowa.edu/arctas/emission.html), fireemissions from FINN-v1 (Wiedinmyer et al., 2011), as well as biogenic emissions based onMEGAN.For both spring and summer model periods, output was downloaded for a wide swath of theglobe (20°N to 80°N latitude (31 cells) and -90°W to -180°W longitude (37 cells)). MOZARToutput contains 56 vertical cells. Output was passed through the program MOZART2CAMx,which reads MOZART files and horizontally and vertically interpolates the data onto a CAMxdomain to generate date-specific initial and/or boundary conditions. Boundary and ICs weregenerated only for the outermost domain, grd01.2.5 Developing Other Model InputsTwo remaining input categories required for CAMx also consist of models, albeit smallerthan WRF, SMOKE, MEGAN and MOZART. These are the photolysis rate models and theland use-leaf area index merging tool. Photolysis rate inputs are critical for the function anyphotochemical model. Final photolysis inputs were processed by two sub-models, O3MAP(which defines the atmospheric ozone column intervals for CAMx) and the TUV (which is aradiative transfer model developed and distributed by NCAR that develops clear-sky photolysisrate inputs (ENVIRON, 2013)). MERGE_LULAI merges independently-developed landuse50and/or leaf area index fields (taken from MEGAN input) with an existing CAMx surface filegenerated by WRFACMx.2.6 CAMx Model SettingsModel settings refers to the various switches and parametrization schemes used in CAMx.Model settings unique to this research are briefly presented below.2.6.1 AdvectionThe Piecewise Parabolic Method (PPM), developed by Colella and Woodward (1984), was se-lected as the horizontal advection solver. This scheme is one of two options available in CAMxand was selected because, despite increasing model runtime, it is considered to have increasedaccuracy over the area-preserving flux-form advection solver developed by Bott (1989).2.6.2 Chemical MechanismThe Carbon Bond (2005 update) (CB05) mechanism was selected for this research. CB05employs 51 species in 156 reactions, and is suitable for computer modelling studies of O3, PM,visibility, acid deposition and air toxics (Yarwood et al., 2005). This is the default option ofthe SMOKE emissions processor and was chosen mainly out of convenience. It is described indetail in Yarwood et al. (2005) (including a list of all reactions and rate constants).2.6.3 Aerosol OptionThe Coarse-Fine (CF) aerosol scheme was selected for this research. The CF scheme divides thePM size distribution into two static modes (coarse and fine). Primary species can be modelledas fine and/or coarse particles, while all secondary species are modelled as fine particles only51(ENVIRON, 2013). This scheme was chosen for simplicity, though aerosols are not the focusof this research.2.6.4 Plume in Grid OptionGrid resolution plays a critical role in resolving emissions advection, yet practical and the-oretical considerations suggest 1000 m as the limit to the resolution of Eulerian air qualitymodels (ENVIRON, 2013). This is an important limitation of the Eulerian-based modellingsystem as it affects their ability to treat point sources. In order to properly model pointsources, sub-models referred to as PiG models were developed (Karamchandani et al., 2011).These sub-models track point source plumes outside of the main Eulerian model (using La-grangian physics and individual plume segments) and only return the plume (also referred toas ’dumping’) to the Eulerian model when:• the size of the puff is approximately the same size as the grid, or• the puff has reached a state of chemical maturity.The Greatly Reduced Execution and Simplified Dynamics (GREASED) PiG option was se-lected for this research. This option is described as ideal for treating the early chemicalevolution of large NOX plumes (ENVIRON, 2013) which matches the emissions profiles ofLNG and marine sources. Also, this option is compatible with aerosol chemistry and the O3source apportionment tools available in CAMx. All stacks greater than 20 m were allotted forplume in grid treatment. While PiG treatment is generally reserved for large sources, a lowthreshold of 0.01 tonnes per day was selected for this research. The result from using such alow threshold was that many NOX sources were immediately diluted and returned into thegrid the following time step.522.6.5 Wet and Dry DepositionWet and dry deposition are important scavenging processes for both aerosols and solublegaseous pollutants. Both wet and dry deposition options were used in this research. TheZhang et al. (2003) scheme was selected for dry deposition as it is described by ENVIRON(2013) as a state-of-the-science method that has been shown to reproduce observed fluxes ofO3 with reasonable accuracy.53Chapter 3Results and Discussion I: Answering Ques-tion 1Recall Question 1 from section 1.5:Can spring and summer photochemical O3 be replicated by a model in a mannerthat is fit for the purpose of investigating worst-case concentrations that may resultfrom proposed industrial emissions in a constrained coastal airshed with complexterrain, such as the Terrace - Kitimat valley?The nature of answering this question lies in the numerous smaller decisions and settings thatcomprise a modelling approach, as outlined in chapter 2. Recreating worst-case photochemicalO3 generally refers to the magnitude of afternoon maxima, however the model would besomewhat inadequate if it could not also broadly recreate time series consistent with the overalldiurnal pattern that O3 exhibits. Also, while O3 is stated explicitly in the above question,some attention must also be given to NO and NO2 because of the important relationshipbetween these gases and O3 as described by reactions 1.1 and 1.2. Therefore, model fitnesstests concerning O3 must also include the fitness of the gases that comprise NOX .54The remainder of this chapter contains an evaluation of model fitness and ultimately conclu-sions regarding the appropriateness of the modelling approach as described in the previouschapter for recreating worst-case spring and summer O3. This chapter is limited to the O3portion of the evaluation while the NO and NO2 portions located in appendices D and E.3.1 Model Fitness Evaluation DescriptionAs ambient monitoring of O3, NO and NO2 is not routine in Kitimat and is limited to two briefperiods which do not coincide with either spring or summer model period, modelled outputwas instead evaluated at the St. Josephs air quality monitoring station in Smithers whereambient observations were available for comparison. Smithers is outside the 1.333 km domainbut is within the 4 km domain; its location relative to Kitimat and Terrace can be seen infigure 2.9. Properties of this station are listed in table 3.1 below.Table 3.1: Location of Smithers St. Josephs air quality monitoring station and correspondingx and y-grid cells in the 4.0 km domains, d03 (WRF) and grd02 (CAMx).Name Lat Lon Elevation Elevation (WRF) x-cell y-cell x-cell y-cellm-asl m-asl (WRF) (WRF) (wrfcamx) (wrfcamx)Smithers St. Josephs 54.783276 -127.177340 497 626 82 73 80 71The model evaluation consisted of three components:1. a qualitative evaluation of modelled and observed time series to check their generalagreement,2. a quantitative evaluation of overall and maximum modelled vs. observed mixing ratiosusing various statistical metrics, and3. a qualitative evaluation of modelled vs. observed scatterplots to visualize their corre-lation and determine whether the model output fall in the 2:1 or 1:2 ratio which iscommonly accepted as a measure of performance for air quality models (Chang andHanna, 2004).Dennis et al. (2010) describe this evaluation as an operational evaluation, where statistical and55graphical analyses are aimed at determining whether model estimates are in agreement with theobservations in an overall sense. The approach to model fitness is similar to those used in thepast (Steyn et al., 2013). The statistical part of the model evaluation was guided by Willmott(1982), Dennis et al. (2010) and Simon et al. (2012). Willmott argued that the best metricsfor model evaluation are difference measures as opposed to correlation measures (Willmott,1982). The strength of these metrics are reinforced by Simon et al. (2012), who listed allstatistical metrics used in evaluating photochemical model performance. Metrics selected forthis evaluation include normalized mean bias (NMB), mean bias (MB), root mean squarederror (RMSE), mean absolute error (MAE), and the correlation coefficient (r). Equations forthese metrics are presented in appendix B. The above were selected as a suite to succinctlyinform model performance based on the following rationale:• NMB allows the bias of various parameters to be compared against each other.• MB is an important element to understanding the performance of any model and wasthe most commonly used metric in Simon et al. (2012).• RMSE is a popular choice in many sciences. It penalizes outliers and gives a goodindication of whether there were many of these in the dataset.• According to Willmott and Matsuura (2005), MAE offers a natural measure of aver-age error as it “avoids the physical artificial exponentiation that is an artifact of thestatistical-mathematical reasoning from which RMSE comes” (Willmott, 1982).• While r is not a difference measure and is not recommended by Willmott (1982), itis nevertheless popular (as documented by Simon et al. (2012)), and adds clarity tomodelled vs. observed scatterplots.Consistent with Dennis et al. (2010), statistical confidence levels are not presented. It shouldalso be noted that there is no pass or fail point for any of these metrics (though Chang andHanna (2004) indicate a NMB < 30% constitutes a ’good’ outcome), rather the model’s fitnessis evaluated based on the cumulative evidence of all metrics.56Dates OmittedNot all dates from the model periods were included in the evaluation. Those removed fromthe analysis are listed in table 3.2, along with the rationale for their exclusion.Table 3.2: Dates removed from model evaluationPollutant Season Date RationaleO3Spring April 14 First day of model run - spinupSpring May 14 Final day of model, run incompleteSummer July 29 First day of model run - spinupSummer August 10 Incomplete observation recordSummer August 15 - 17 Forest fireSummer August 18 Final day of model run incompleteNO and NO2Spring April 14 First day of model run - spinupSpring May 2 - 3 Incomplete observation recordSpring May 14 Final day of model run incompleteSummer July 29 First day of model run - spinupSummer August 15 - 17 Forest fireSummer August 18 Final day of model, run incomplete3.2 Evaluation Results3.2.1 OzoneAs can be seen in figures 3.1 (a) and 3.1 (b), modelled O3 and observed mixing ratios showgeneral agreement during daytime hours. Maximum modelled O3 is similar to maximummeasured O3 on most days in the spring, while the summer performance is slightly poorer.Overnight however, the model was unable to capture low O3 values. The observed meandifference between maximum and minimum O3 mixing ratios were 36.8 ± 1.6 ppb in thespring and 28.2 ± 1.8 ppb in the summer, while the modelled mean difference was only 6.9± 0.6 ppb in the spring and 8.8 ± 1.3 ppb in the summer (note: ± values represent standarderror).57Figure 3.1: Observed and modelled O3 time series at Smithers for: a) spring and b) summermodel periods. Black lines are observed and red lines are modelled Control Case.Evaluation statistics are listed in tables 3.3 and 3.4. Table 3.3 refers to all model-observedhours, while table 3.4 refers to a 12 hour period from 10:00 to 21:00 (the last hour refersto data from 21:00 to 22:00). Considering all paired hours (top line of table 3.3), statisticalmetrics indicate a clear high bias in both spring and summer, with a MAE of 15.4 and 12.1 ppbrespectively. Considering maximum modelled vs. maximum observed values (middle line oftable 3.3), statistical metrics are much better. A comparison of the maximum daily modelledmixing ratio and the value of its paired observed value (bottom line of 3.3) indicate CAMx isnot able to predict the correct hour of maximum observed O3.58Table 3.3: O3 model evaluation statistics for the Control case - (all hours)Spring SummerDescription n NMB MB RMSE MAE r n NMB MB RMSE MAE r(%) (ppb) (ppb) (ppb) (%) (ppb) (ppb) (ppb)Paired mixing ratios 662 51.672 13.989 19.921 15.396 0.108 344 63.834 11.072 15.706 12.095 0.189Max model / max obs 29 2.543 1.097 4.857 3.862 0.593 15 5.203 1.622 6.169 4.885 0.235Max model w/ paired obs 29 63.027 17.104 22.514 17.811 0.209 15 46.975 10.482 16.174 11.706 -0.289Statistical metrics improved when data from only the 12 hour period from 10:00 to 21:00 wereused. When considering all paired modelled vs. observed values spring and summer biasesdecreased from 51.7 to 12.8% and 63.8 to 18.9% respectively, RMSE decreased from 19.9 and15.7 ppb to 10.2 and 9.2 ppb (spring and summer respectively) while MAE decreased from15.4 and 12.1 ppb to 7.4 and 6.6 ppb (spring and summer respectively). These biases are wellbelow the 30% value presented by Chang and Hanna (2004) as ’good’. The metrics for dailymaximum modelled vs. maximum observed mixing ratios also improved, with a MB of lessthan 1 ppb and a reasonable r value (0.65) in the spring.Table 3.4: O3 model evaluation statistics for the Control case - (10:00 to 21:00)Spring SummerDescription n NMB MB RMSE MAE r n NMB MB RMSE MAE r(%) (ppb) (ppb) (ppb) (%) (ppb) (ppb) (ppb)Paired mixing ratios 333 12.791 4.693 10.219 7.415 0.265 174 18.923 4.665 9.214 6.568 0.371Max model / max obs 29 1.597 0.683 4.787 3.772 0.649 15 2.332 0.727 4.871 4.071 0.523Max model w/ paired obs 29 30.441 10.145 15.939 11.315 0.270 15 20.774 5.487 9.216 7.738 0.342Figure 3.2 is a scatter plot of all 10:00 to 21:00 spring (red dots) and summer (blue dots) pairedobserved and modelled maximum mixing ratios, along with observed maximum and modelledvalues for the spring (red diamonds) and summer (blue circles). This figure demonstrates thestrength of the Control case output, which is that maximum modelled O3 mixing ratios aresimilar to maximum observed. Given the intent of research question 1, agreement betweenmaximum modelled and observed values is crucial.On the other hand, better agreement between modelled and observed O3 over all hours wouldimprove model fitness and better demonstrate the appropriateness of the model approachfor answering Question 1. Among all the reasons for the inability of CAMx to successfullyreproduce overnight O3 mixing ratios, two were investigated:59Figure 3.2: Observed vs. modelled O3 scatterplot at Smithers for spring and summer modelperiods. Small red diamonds small blue dots are daytime (10:00 - 21:00) hourly paired valuesfor spring and summer Control Case model periods respectively. Large red diamonds and largeblue circles are maximum daytime modelled vs. maximum daytime observed mixing ratios forspring and summer Control Case model periods respectively.1. emissions were improperly characterized, and2. emissions diffused out of the surface layer to higher layers in the model overnight, re-ducing NOX available for titration.The complete evaluation of NO and NO2 model output assisted in addressing the first pointand an examination of output through vertical layers addressed the second point. Otherpossible causes of the lack of agreement between overnight modelled and observed mixingratios, include: the failure of WRF to adequately characterize some meteorological parametersand potentially CAMx’s inability to adequately characterize nocturnal chemistry in rural areas.These were not investigated, as the NO and NO2 time series, as presented in appendix D,clearly indicate that emissions were a large contributor to the problem. Vertical NO and NO2plots presented in appendix D for Smithers also demonstrate that diffusion did not contributeto the poor results overnight.603.2.2 Fitness of Control CaseAs seen from the time series plots, O3 scatter plot and statistical metrics, the control caseemissions as described in section 2.3.1, captured afternoon peak O3 values reasonably but failedto recreate the diurnal O3 pattern that was measured in Smithers. Also, as shown in appendixD, the model failed to replicate NO and NO2 mixing ratios at all hours of the day. Indeed, itappeared as though CAMx acted independent of the anthropogenic NOX emission inventory.Given that the desired outcome of this research was to investigate not only peak O3 but alsothe complete spatiotemporal O3 distribution, the model’s original emissions inventory causedCAMx to be unfit for research purposes. A series of sensitivity model runs was conducted todetermine a more appropriate baseline emission inventory, which is the subject of the nextsection.3.3 CAMx Sensitivity to Mobile and Point Source EmissionsPerturbations3.3.1 Sensitivity Case StudiesA series of sensitivity tests were conducted as part of the model evaluation to gauge whetherimprovements to model fitness could be attained for O3, NO and NO2 by perturbing NOXemissions. Sensitivity cases can be considered a form of top-down emission inventory develop-ment, where the final concentrations are known and the inventory is built to match what wasobserved. Five scenarios were selected for the sensitivity study, and are listed and described intable 3.5. Emissions for a single cell in Smithers are listed, along with emissions for Smithersand the surrounding cells. The final column lists valley-wide NOX emissions in tonnes perday. It should be noted that only NOX emissions were perturbed and that no other CAC wasmodified.Three different approaches to emissions perturbations were investigated:61Table 3.5: Sensitivity test names and totalled NOX emissions in SmithersScenario One Cell Emissions One Cell All Smithers Emissions Approx Smithers TotalDescription Name Mobile Point Bio (moles NO/hr) Mobile Point Bio (moles NO/hr) (moles NO2/hr) (tonne NOX/day)Base Base 38.1 103.1 1.0 141.2 343.1 206.2 10.0 559.3 62.1 0.47110x Mobile base Point 10xm 381.2 103.1 1.0 484.3 3430.8 206.2 10.0 3647.0 405.2 3.074Base Mobile 20x Point 20xp 38.1 2026.7 1.0 2065.8 343.1 4053.4 10.0 4406.5 489.6 3.71410x Mobile 20x Point 10xm20xp 381.2 2026.7 1.0 2418.9 3430.8 4053.4 10.0 7494.2 832.7 6.31750x Mobile base Point 50xm 1905.7 103.1 1.0 2019.8 17154.0 206.2 10.0 17370.2 1930.0 14.642100x Mobile 20x Point 100xm20xp 3812.0 2026.7 1.0 5839.7 34308.0 4053.4 10.0 38371.4 4263.5 32.3441. perturbing mobile, time varying emissions,2. perturbing point, time invariant emissions, and3. perturbing both mobile and point sources.Three of the five cases (10xm, 20xp and 10xm20xp) represent modest emissions perturbations(relative to the original emissions inventory), while two (50xm and 100xm20xp) representextreme cases that are unrealistic in magnitude. These emission perturbations were chosenfor a variety of reasons. Time invariant point source emissions were perturbed by 20 timesto investigate the model’s ability to titrate O3 overnight. Recall that the point sources listedby the NPRI for Smithers had no stack information and thus were treated as area sources byCAMx. Thus these increases in NOX can be considered to account for any and all sourceswhich are constant throughout the day (such as industrial emissions unaccounted for by theNPRI). Mobile emissions were perturbed by a variety of magnitudes (10, 50 and 100 times)to investigate the model’s response to emissions occurring throughout the length of the entireBV.The diurnal weighting of the time-varying emissions is presented in figure 3.3. From this it canbe seen that mobile emissions are greatest in the morning and evening hours, with minimalemissions overnight.O3 results from the sensitivity tests are presented in the following subsections while the resultsfor NO and NO2 are presented in appendix E. Results for O3 demonstrate improved modelfitness with increasing emissions. CAMx had difficulty matching the timing of peak NO valueswhich led to poor statistical outcomes for that pollutant and NO2 even though the magnitudesof peak values for improved noticeably.62Figure 3.3: Hourly distribution of daily point and mobile emissions. Red line is the hourlydistribution of point source emissions while the blue line is the hourly distribution of mobilesource emissions.3.3.2 Ozone Sensitivity ResultsAs can be seen from the figures 3.4 (a) and 3.4 (b), nighttime titration was best emulated bythe most extreme case, 100xm20xp. The 50xm case also led to overnight titration though itsperformance was inconsistent (not surprising as this case had reduced emissions overnight).Emissions from the 20xp and 10xm20xp cases also led to some reduced overnight O3 whilethe 10xm case made little difference to overnight O3. Daytime O3 maxima were increased inall cases, at times bringing the maxima closer to the observed maxima and at times doing theopposite.Based on the statistical metrics listed in tables 3.6 and 3.7 it is evident that for O3, modelperformance improved with increasing NOX emissions. Comparing the Base case to the100xm20xp case, in the spring model bias decreased from 51.7% to 14.2% (all hours, table3.6) and from 12.8% to 2.3% (daytime hours, table 3.7). In the summer, model bias decreasedform approximately 63.8% to 14.8% (all hours, table 3.6) and from 18.9% to 4.0% (daytimehours, table 3.7). The r also increased with each case of increasing emissions, from 0.11 to 0.62(spring Base to 100xm20xp, table 3.6) and from 0.19 to 0.71 (summer Base to 100xm20xp,63Figure 3.4: Observed and modelled O3 time series at Smithers for: a) spring and b) summermodel periods. Black lines are observed and red lines are modelled Control Case. Green,blue, purple, orange and light blue lines are the modelled sensitivity test cases: 10xm, 20xp,10xm20xp, 50xm and 100xm20xp respectively, as defined in table 3.5. 64table 3.6). Considering all cases, springtime bias metrics were lower than summertime biasmetrics, while error metrics were lower in the summer than the spring. Correlation was alsobetter in the summer period (all cases). For O3, increased emissions led to improvements inmodel fitness for all cases except for the 10xm case, which did little to reduce overnight O3and led to poorer performance statistically due to increased daytime O3 mixing ratios.Table 3.6: O3 model evaluation statistics for the sensitivity tests - (all hours)Spring SummerDescription n NMB MB RMSE MAE r n NMB MB RMSE MAE r(%) (ppb) (ppb) (ppb) (%) (ppb) (ppb) (ppb)Base 662 51.672 13.989 19.921 15.396 0.108 344 63.834 11.072 15.706 12.095 0.18910x Mobile base Point 662 51.683 13.992 19.404 15.133 0.270 344 69.323 12.024 16.012 12.611 0.343Base Mobile 20x Point 663 48.820 13.230 18.595 14.537 0.352 344 52.369 9.084 13.136 10.250 0.46710x Mobile 20x Point 662 42.245 11.437 17.139 13.141 0.390 344 60.170 10.437 13.897 11.115 0.52250x Mobile base Point 664 45.638 12.370 17.295 13.775 0.509 344 58.364 10.123 13.480 10.846 0.599100x Mobile 20x Point 664 14.170 3.841 12.027 9.353 0.624 344 14.822 2.571 8.656 6.527 0.706Table 3.7: O3 model evaluation statistics for the sensitivity tests - (10:00 to 21:00)Spring SummerDescription n NMB MB RMSE MAE r n NMB MB RMSE MAE r(%) (ppb) (ppb) (ppb) (%) (ppb) (ppb) (ppb)Base 333 12.791 4.693 10.219 7.415 0.265 174 18.923 4.665 9.214 6.568 0.37110x Mobile base Point 334 15.467 5.675 10.612 7.862 0.313 174 27.440 6.765 10.524 7.811 0.431Base Mobile 20x Point 335 14.877 5.462 10.153 7.779 0.408 174 14.382 3.546 7.408 5.664 0.60210x Mobile 20x Point 333 10.822 3.970 9.477 7.062 0.385 174 23.278 5.739 8.879 6.24 0.58750x Mobile base Point 335 19.141 7.023 11.447 8.938 0.428 174 29.945 7.383 10.501 8.390 0.601100x Mobile 20x Point 335 2.328 0.854 9.378 7.225 0.503 174 4.011 0.989 6.825 5.339 0.749Finally, figure 3.5 shows an O3 scatter plot of the 20xp case. This case was chosen to demon-strate how, when emissions are increased in a modest fashion, the strength of the original modelresults (reasonable recreation of maximum mixing ratios) is not compromised. Included areall 10:00 to 21:00 spring (red dots) and summer (blue dots) paired observed - modelled mixingratios along with observed maximum and modelled values for the spring (red diamonds) andsummer (blue circles). Considering how overnight O3 titration is increased with this case, theaddition of a constant low-level NOX area source is seen as improving model fitness.These additional emissions could conceivably represent two potential NOX sources in Smithersthat were unaccounted in the original emissions inventory:65Figure 3.5: Observed vs. modelled O3 scatterplot at Smithers for spring and summer modelperiods. Small red diamonds small blue dots are daytime (10:00 - 21:00) hourly paired valuesfor spring and summer 20xp sensitivity case model periods respectively. Large red diamondsand large blue circles are maximum daytime modelled vs. maximum daytime observed mixingratios for spring and summer 20xp sensitivity case model periods respectively.• Residential and commercial heating.Residential and commercial heating in the form of natural gas (or wood) combustionis a source of NOX emissions not included in the emissions inventory. Particularly forthe spring period when overnight temperatures in Smithers often fall below 0 °C, it isnot unreasonable to assume that these emissions occur overnight and contribute to O3titration.• Other Industrial sources unaccounted for in the NPRI.The NPRI contained one emission source for each Smithers-based industry, however inreality these facilities were comprised of numerous emission sources. For example, actualemissions at the sawmill included natural gas-fired kilns (Pierce, 2015), yet these werenot identified in the NPRI. Also, the particleboard plant contained two wood-fired dryingstacks (which emit NOX) however the NPRI only listed one entry for that facility. Both66of these facilities were located upwind of Smithers at night (when the planetary boundarylayer is stable). Discrepancies in emissions from industrial facilities could explain whyincreasing low-level emissions increased model fitness, as these sources should have beenincluded in the NPRI’s original emission inventory.The results of the NO and NO2 sensitivity tests can be found in appendix E, along withscatterplots for each sensitivity test. In this appendix it is shown that NO and NO2 results alsoimproved with increasing emissions, however statistical results for NO2 were most improvedfor the 20xp and 10xm cases, as excess NO emissions let to overproduction of NO2 for the50xm and 100xm20xp cases.3.4 Question 1 Summary and ConclusionThe following is a summary of the salient outcomes of the five sensitivity model runs.• The emission inventory for Smithers, taken from the NPRI did not adequately captureNOX emissions in Smithers. This was evident from the Control case model run where O3output had little diurnal variability and from NO output which, as shown in appendixD, rarely exceeded 2 ppb.• CAMx was capable of titrating O3 overnight given additional NO emissions. As shownin appendix E, overnight O3 mixing ratios decreased with increasing NO emissions.• As shown in appendix E, when NO emissions were very high, overnight NO2 productionwas elevated and became unrealistically high.• Modest increases in NOX emissions from the 10xm and 20xp scenarios led to improvedstatistical metrics for all gases (O3, NO and NO2) when compared with the Base case.• With the exception of April 21, none of the modest increase scenarios resulted in greaterthan 25% elevated O3 when compared with the Base case, suggesting that daytimemaxima for these cases were similar with the Base case.67• Of the modest-increase scenarios, the 20xp scenario was better able to titrate overnightO3 then the 10xm case.Based on the above, the following conclusion can be drawn:• A top-down emission inventory approach can be a reasonable method for estimatingemissions for the BV, though a well-developed spatial and temporal emissions inventorywould have been ideal. Care must be taken to ensure a balance between elevated emis-sions (which improves O3 and NO) and elevated titration (which produces too muchNO2).Considering the above, an affirmative answer can be assigned to Question 1. With the inclusionof an updated emissions inventory, photochemical O3 was replicated by a model in a mannerfit for the purpose of investigating worst-case concentrations. The emission inventory for bothTest and Control cases in Kitimat was subsequently amended to incorporate the results of themodel evaluation. In addition to all proposed sources in Kitimat a low-level emission sourcewas added to the inventory in order to facilitate overnight O3 titration. This approach doesnot change maximum daytime O3 (which forms downwind of NOX emissions (Sillman, 2012))and creates more realistic results overnight, as seen in the above sensitivity tests. In orderto prevent over-titration of O3 the new source was minimal, approximately 50% of the 20xpcase (˜1.5 tonnes of NOX per day). These emissions were uniformly distributed throughoutthe urban community of Kitimat and had no diurnal variability. This change to the emissioninventory was a top-down approach to account for all unaccounted sources in the originalemissions inventory. This additional source, when contrasted with the all proposed NOXsources in Kitimat, is minor in magnitude, as seen from the small purple addition to theNOX column in figure 3.6. Also, an updated image showing emission sources near Kitimat ispresented in figure 3.7.68Figure 3.6: Updated daily NOX emissions in the Terrace - Kitimat valley for the Test case.Other emissions remain unchanged.69Figure 3.7: Updated emission locations in Kitimat for the Test case. identifies locationsof aluminum smelter sources, identifies locations of liquefied natural gas sources andidentifies locations of marine sources. ￿ identifies the locations of new low-level area NOXsources. Parallel lines and the solid black line indicate locations of mobile sources (roads andrail sources respectively). Note that local roads are identified with thin brown lines thoughthese are not associated with mobile emissions.70Chapter 4Results and Discussion II: AnsweringQuestion 2Question 2 has two parts, both of which are useful in deepening our understanding of worst-caseO3 production in the TKVA. The exact questions is: should all proposed industrial facilitiesin the TKVA be constructed,(a) Where are the locations of the O3 maxima, and what are the worst-case concentrationsthat could occur there?(b) What meteorological factors contribute to enhancing or reducing O3 concentrations inthe Terrace - Kitimat valley airshed?Because of averaging times and statistical metrics a number of meanings could be assigned tothe words worst-case. In this research worst-case is thought of as hourly maxima, though asO3 research evolves in the TKVA other forms of worst-case may be considered. An importantthreshold for this analysis is the one-hour mixing ratio of 82 ppb, the level at which theprovincial government issues O3-based air quality advisories (British Columbia Ministry ofEnvironment, 2014).71This section follows the following format: section 4.1 looks at the results of the Test casefor the spring period. Time series for the three TKVA communities are presented. Spatialand temporal maxima and minima are explored, along with an analysis of the meteorologicalconditions contributing to the greatest change in O3 mixing ratios. Section 4.2 is similar innature to 4.1, though instead of presenting time series at all communities, time series are pre-sented for Kitimat and the locations of the spatial maxima. An analysis of the meteorologicalconditions contributing to the greatest change in O3 mixing ratios is included, as well as anexample of how varied initial conditions can affect model results.4.1 Springtime4.1.1 Time Series: Communities in the Terrace - Kitimat Valley AirshedCommunity time series plots tell one part of the spring model period’s story. The springmodel period was selected to investigate whether increased NOX emissions would exacerbateO3 formation during the springtime peak which is measured across mid and high latitudes inthe northern hemisphere (Monks, 2000). Test Case results for the communities of the TKVAindicated that in the spring elevated NOX emissions often led to the opposite, and contributedto broad O3 titration and NO2 production. This scavenging was most prevalent in the valleybottom overnight and was shown to, at times, decrease O3 mixing ratios by 50%.Figure 4.1 contains spring period O3 time series for Terrace, Kitimat and Hartley Bay. Includedin these plots are planetary boundary layer height (PBLH) and temperature time series. Theshaded areas indicate a south wind, while the non-shaded areas indicate a north wind. Ascan be seen from these plots, Test Case emissions made little difference to O3 mixing ratiosin Hartley Bay (fig 4.1 (a)). Closer to the emission sources in the community of Kitimat,excess (fig 4.1 (b)) O3 was also rarely produced. Indeed during the model period O3 was oftenscavenged, particularly when the wind was from the south. Reductions in O3 came at theexpense of NO2 production (not shown). O3 titration was amplified from April 15th - 18th72when there were horizontally recirculation winds over the course of several days. As shown bythe alternating shaded and non-shaded backgrounds of the time series, one can imagine the’sloshing’ of air up and down the valley with this repeated wind pattern. Lastly, O3 mixingratios in Terrace were largely unchanged by the addition of new industrial emissions, exceptfor April 16th and 17th, when there was some overnight O3 scavenging by the NOX plume thatstretched the entire length of the TKV (fig 4.1 (a)).4.1.2 Locations of Minima and MaximaSpatially the results are more nuanced. Overnight O3 titration was often prevalent at lowelevations along the entire length of the TKV when winds were southerly. This is explainedby the large-scale NOX sources and the shallow boundary layer which trapped these plumeswhen mixing was at a minimum. As an example, figure 4.2 (a) displays an 80 km zone ofreduced O3, while figure 4.2 (b) relates these mixing ratios as a percent change from theControl case. Because the marine emission sources stretch the length of the Douglas channel,O3 is scavenged along the entire length of the CAMx grd03 domain (with the exception of thefar-north section).73Figure 4.1: Planetary boundary layer depth, O3 and temperature time series for the springmodel period. Model output at: a) Terrace, b) Kitimat and c) Hartley Bay. Black lines arethe Control case and red lines are the Test case. Purple shading indicates southerly wind asoutput from the WRF in Kitimat.74Figure 4.2: CAMx spatial output for the Terrace - Kitimat valley airshed on April 17th, 2010 at 00:00 PST showing: a) Test Casemodelled O3 mixing ratios, and b) percent change from the Control case. The green line outlines the Terrace - Kitimat valleyarished. ￿, ￿ and ￿ identify the communities Terrace, Kitimat and Hartley Bay respectively, while and identify hotspotlocations Lakelse Lake and Miskatla Inlet. Parallel lines identify the locations of roads.75On the other hand, daytime production of O3 also occurred, albeit in modest quantities. Thehours with the greatest O3 production were 12:00 - 16:00 PST, with 14:00 PST being the hourwith greatest O3 differences between Test and Control cases. Locations of this productionvaried with wind direction but were typically along valley walls, downwind of emissions andoutside the fresh NOX plume. In absolute terms, the change in ambient O3 mixing ratios waslimited in magnitude to a maximum of 5 ppb (and often less) though spatially these increasesoccurred over broad sections of the domain. As the spring O3 peak diminished over the monthof April and into May these modest increases in O3 led to a greater overall percent changein O3. Figure 4.3 (a) contains an image showing O3 mixing ratios at 16:00 PST on May 4th,while figure 4.3 (b) shows the percent difference between these mixing ratios and the Controlcase.76Figure 4.3: CAMx spatial output for the Terrace - Kitimat valley airshed on May 4th, 2010 at 16:00 PST showing: a) Test Casemodelled O3 mixing ratios, and b) percent change from the Control case. The green line outlines the Terrace - Kitimat valleyarished. ￿, ￿ and ￿ identify the communities Terrace, Kitimat and Hartley Bay respectively, while and identify hotspotlocations Lakelse Lake and Miskatla Inlet. Parallel lines identify the locations of roads.774.1.3 Ingredients for Reduced OzoneThe meteorological factors which contributed to the greatest change in O3 over the springperiod were the coupling of low PBLHs with diurnal horizontal recirculation, particularly onthe longer nights early in April. These variables, which contributed to O3 production in thesummer (with the exception of night length, see section 4.2), led to long-lasting titration of O3along the TKV bottom in the spring. Figure 4.4 displays spring PBLH at four locations in theTKV along with generalized mountain peak elevations (as determined from contour plots).Figure 4.4: Planetary boundary layer depth time series for the spring model period in theTerrace - Kitimat valley airshed. Model output at: Terrace (yellow lines), CYXT (orangelines), Kitimat (blue lines) and Hartley Bay (green lines). Static lines indicate nearby mountainpeak elevations, coloured for the respective communities.As can be seen from this image in conjunction with the time series in figure 4.1, the dayswith low PBLH corresponded to days with increased differences between the Control andTest case O3 mixing ratios. This was amplified on days with low PBLH as well as horizontalrecirculation, exemplified for the period of April 15th - 19th. As can be seen by the shadedsections in figure 4.1, there was a 5-day diurnal onshore-offshore wind pattern which advectedpollutants up and down the valley. During this time the low PBLH acted like a lid and limiteddispersion vertically. However, despite contributing to an environment for increased O3 as78occurred during the summer (see section 4.2), this created an environment for long-lasting O3titration. These dates in the mode period have the longest nights and fewest hours of daylight,another possible contributor. (In fact, these are the shortest days of either model period.) Itis also possible that despite ever increasing NO2 concentrations, the valley lacked sufficientVOC emissions to lead to meaningful O3 production as occurs in reactions 1.5 through 1.9. Asseen in figure 2.14 (a), isoprene emissions were minimal during this time (owing to the absenceof leaves on shrubs and deciduous trees) and while terpene emissions were dominant, the ratecoefficient of ↵-Pinene (an example terpene) against OH is 50% of isoprene’s (Derwent, 1999).Thus, in the absence of VOC emissions, those factors which create O3 can also lead to reducedO3.4.2 Summertime4.2.1 Time Series: Hotspots in the Terrace - Kitimat Valley AirshedThe summer period resulted in increased O3 production in certain locations of the TKVAunder certain meteorological conditions, often in the 12 - 20 ppb range. A spatial examinationof the results shows two two prevalent ’hotspots’ (i.e.: areas with the greatest increase overthe Control case), one to the south of the emission sources and the other to the north.To the south, the location with the greatest O3 increase was Miskatla Inlet, a small bay 30 kmsouth of Kitimat along the Douglas channel. Miskatla Inlet was not a location with high VOCemissions; instead this location seemed to trap recirculating air which had advected southward(from Kitimat) earlier in the day and was returning north under an onshore wind. As seen infigure 4.5 (c), Test Case O3 at this location was at times 50% higher than the control case.To Kitimat’s north the location of the greatest O3 increase was Lakelse Lake, 13 km south ofTerrace. As seen in figure 4.5 (a), the difference in Control and Test Case O3 at this locationreached over 55% on one instance, and is the subject of further analysis in section 4.2.3. Asseen in figure 4.5 (b), O3 increases did not occur (in any meaningful way) in the community of79Kitimat itself because, as discussed earlier, this pollutant manifests downwind of the emissionsof its precursors (Sillman, 2012).4.2.2 Locations of Minima and MaximaWhile the locations of greatest increase were within the 1.333 km domain, the spatial ex-tent of O3 production extended beyond its boundaries on dates which were conducive to O3development. Nevertheless only results for the 1.333 km domain are presented. O3 produc-tion occurred throughout the Douglas channel as well as the lower TKV, depending on winddirection. Under a north wind, air masses advected south through the Douglas channel. Un-der a pattern of horizontally recirculating wind (as exhibited during the first 5 days of thesummer period), as these air masses returned northward, pollutants were funnelled into andtrapped in Miskatla Inlet, explaining why this location to the south of Kitimat is the locationof the southern maximum O3. Figures 4.6 (a) and 4.6 (b) show the accumulation of O3 there(coloured with a blue circle) as well as the percent increase from the Test case over the Controlcase.80Figure 4.5: Planetary boundary layer depth, O3 and temperature time series for the summermodel period. Model output at: a) Lakelse Lake, b) Kitimat and c) Miskatla Inlet. Blacklines are the Control case and red lines are the Test case. Purple shading indicates southerlywind as output from WRF in Kitimat.81Figure 4.6: CAMx spatial output for the Terrace - Kitimat valley airshed on August 4th, 2010 at 15:00 PST showing: a) TestCase modelled O3 mixing ratios, and b) percent change from the Control case. The green line outlines the Terrace - Kitimat valleyarished. ￿, ￿ and ￿ identify the communities Terrace, Kitimat and Hartley Bay respectively, while and identify hotspotlocations Lakelse Lake and Miskatla Inlet. Parallel lines identify the locations of roads.82Under a south wind, pollutants advected northward towards Terrace. The location of thenorthern maximum was adjacent to Lakelse Lake, approximately 13 km south of Terrace andclose to the CYXT airport. Areas surrounding Lakelse Lake have land cover conducive toisoprene emission, as illustrated in figure 2.15. These findings are in general agreement withthose of Castell et al. (2010), who found O3 maxima between 40 and 140 km away from thehypothetical emission source, and who also found that increases were greatest in locationswith high bVOC and low NOX emissions. Lakelse Lake is not a community itself though boththe eastern and western shores are inhabited.Figures 4.7 (a) and 4.7 (b) show the accumulation of O3 at this location (coloured with apurple circle), as well as the percent increase from the Test case over the Control case. Moreinformation on this event is discussed in later sections. It should be stated that at no time didO3 mixing ratios exceed the provincial one-hour objective of 82 ppb. The eight-hour objectiveof 63 ppb was also not exceeded (ignoring the 99th percentile metric shown in table 1.1) eitherat any particular location or following the plume as it advected down then up the valley.83Figure 4.7: CAMx spatial output for the Terrace - Kitimat valley airshed on August 15th, 2010 at 16:00 PST showing: a) TestCase modelled O3 mixing ratios, and b) percent change from the Control case. The green line outlines the Terrace - Kitimat valleyarished. ￿, ￿ and ￿ identify the communities Terrace, Kitimat and Hartley Bay respectively, while and identify hotspotlocations Lakelse Lake and Miskatla Inlet. Parallel lines identify the locations of roads.844.2.3 Ingredients for Elevated OzoneThis section provides an in-depth examination of the elevated O3 as modelled from August 14thto August 15th. During this time it is proposed that the cause of the elevated concentrationswere warm temperatures (which led to high biogenic emissions), low mixing heights (PBLHlower than the height of the surrounding mountain peaks), differential land and water surfaceheating and a two-day period of horizontal recirculation.Warm Temperatures (High Biogenic Emissions)As can be seen from the WRF evaluation and figure C.1 (b) in appendix C, maximum mod-elled temperature increased daily from the afternoon of August 9th through the afternoon ofAugust 14th (with the exception of August 12th which was modelled as slightly warmer thanthe 13th), exceeding 29 °C in Kitimat and also in the TKV. This increase in temperature cor-responded with an increase in overall biogenic emissions in the TKV, as shown in figure 2.14(b). Interestingly, valley-wide bVOC emissions on August 14th were 100 tonnes higher thanAugust 13th even though the maximum temperature only increased by 5 °C between thosetwo days. A large increase in biogenic emissions also occurred from August 3rd to 4th whenmaximum modelled temperature exceeded 29 °C, suggesting that 29 degrees is a good triggerfor high biogenic emissions in the summertime. Figure 4.8 (a) is a linear scatterplot of dailybiogenic emissions vs. maximum modelled temperature for both spring and summer periods.When considered on a seasonal basis, the relationship between maximum temperature andbVOC emissions appear linear, with the summer maximum temperature explaining over 93%of the variability in VOC emissions. However, when both seasons are considered together therelationship appears exponential. Figure 4.8 (b) shows the same points as 4.8 (a) though witha logarithmic-scaled ordinate. From this plot the exponential dependence of bVOC emissionswith temperature is evident, where temperature explains over 95% of the variability in bVOCemissions. This curve matches general results found by Guenther et al. (1993). (Note: above40 °C emissions do begin to decrease)85Figure 4.8: Maximum daily temperature vs. bVOC emissions scatterplots in the Terrace -Kitimat valley airshed for spring and summer model periods presented on: a) linear axesand b) semi-log axes with best fit lines. Blue dots are spring model period and red dots aresummer model period. Individual linear fits for both periods are shown in a) and a combinedexponential fit is shown in b). r2 values are: 0.66 for the spring period, 0.93 for the summerperiod and 0.95 for both periods combined.Low Mixing Heights (PBLH Lower than Surrounding Mountains)Figure 4.9 shows modelled PBLH at four locations in the TKVA for the entire summer period,along with generalized mountain peak elevation (as determined from contour plots). As can beseen in this figure, summertime mixing heights rarely exceeded the height of the surroundingmountain peaks. This led to the accumulation of pollutants, as the only effective outletsduring these times were the north or south ends of the valley. Unfortunately, as the nearestradiosonde is hundreds of kilometres away it was not possible to evaluate modelled PBLH withactual observations.As can be seen, in figure 4.9, PBLHs were low on August 14th and 15th, with maxima of 550m in Kitimat, 1000 m at CYXT and 1200 m in Terrace on the 15th, the day with the highestO3 mixing ratios.86Figure 4.9: Planetary boundary layer depth time series for the summer model period in theTerrace - Kitimat valley airshed. Model output at: Terrace (yellow lines), CYXT (orangelines), Kitimat (blue lines) and Hartley Bay (green lines). Static lines indicate nearby mountainpeak elevations, coloured for the respective communities.Differential Surface Heating and Horizontal RecirculationOver the course of August 14th and 15th differential surface heating over land and wateroccurred along with a repeated pattern of onshore - offshore wind. Differential surface heatingwas prevalent throughout most of the model period, while the repeated onshore - offshore windoccurred in late July and early August, and then again on August 10th, 14th and 15th. While ithas not been established that a sea breeze circulation was present over these days, differentialheating and onshore - offshore wind are two of the main ingredients of a sea breeze, which areknown to coincide with and indeed intensify air pollution episodes (Steyn, 2003). Figure 4.10shows the difference between 2 m land and 2 m sea temperature for the entire summer modelperiod, as output by WRF. The shaded periods indicate a south wind while the non-shadedperiods indicate a north wind. Interestingly, the difference in land and sea temperature wasgreater in early July than in mid August, however the maximum daily temperatures were onlyin the mid 20’s (refer to figure 4.5). This underscores the importance of the presence of allO3-generating ingredients.In the absence of formally declaring this as an occasion where a sea breeze occurred, it is rea-87Figure 4.10: 2m temperature difference between Kitimat and the Douglas channel as modelledby WRF for the summer model period. Purple shading indicates southerly wind as outputfrom WRF in Kitimat.sonable to label this period as one with horizontal recirculation based on the onshore - offshorewind that can be seen in figures 4.11 and 4.12. This pattern recirculates pollutants throughoutthe valley, exacerbating air pollution problems. Figures 4.11 and 4.12 show streamlines (red)in the 1.333 km domain on August 14th (14:00 through 20:00) and 15th (10:00 through 16:00)respectively, at two hour intervals. These figures also show wind barb cross sections along theTKV at 54.1 and 54.4°(the approximate locations of Kitimat and the CYXT airport).These hours were selected to highlight the shift from offshore to onshore wind. Using thecombination of streamlines (which show the general anabatic and katabatic flow within thevalley and other tributary valleys) along with the wind barbs (which provide additional infor-mation about wind speed and direction within the TKV itself), one can see the timing of thetransition (16:00 on August 14th and 12:00 on August 15th) and appreciate how on those daysair would have ’sloshed’ up and down the valley. Figure 4.10 suggests that onshore breezescommence around the time where Tland - Tsea is at a maximum and persist until Tland - Tseais approximately 0 (i.e.: the land has cooled to below the temperature of the nearby ocean).88Figure 4.11: Modelled wind streamlines and wind barb cross sections in the Terrace - Kitimatvalley on August 14th at: a) 14:00, b) 16:00, c) 18:00, and d) 20:00 PST. Red lines arestreamlines and white lines are wind barbs. Grey triangles along the wind barb cross sectionsindicate the midpoint of the cross section.89Figure 4.12: Modelled wind streamlines and wind barb cross sections in the Terrace - Kitimatvalley on August 15th at: a) 10:00, b) 12:00, c) 14:00, and d) 16:00 PST. Red lines arestreamlines and white lines are wind barbs. Grey triangles along the wind barb cross sectionsindicate the midpoint of the cross section.90Spatiotemporal O3 EvolutionBased on the above it is proposed that warm temperatures, low mixing heights, differentialland and water surface heating and a two-day period of horizontal recirculation in combinationprovide the necessary ingredients for elevated O3 mixing ratios in the TKVA. Consideringthat O3 was not produced during periods where one ingredient was missing (e.g.: August4th - August 5th which was missing horizontal recirculation or April 15th - 19th which wasmissing warm temperatures), the importance of all these factors acting in concert becomesclear. Figures 4.13 and 4.14 show the spatiotemporal evolution of O3 on August 14th andAugust 15th respectively.91Figure 4.13: CAMx spatial O3 output in the Terrace - Kitimat valley airshed for the Test case on August 14th, 2010 at: a) 08:00, b)10:00, c) 12:00, d) 14:00, e) 16:00, f) 18:00, g) 20:00 and h) 22:00 PST. The green line outlines the Terrace - Kitimat valley arished.￿, ￿ and ￿ identify the communities Terrace, Kitimat and Hartley Bay respectively, while and identify hotspot locationsLakelse Lake and Miskatla Inlet. Parallel lines identify the locations of roads.92Figure 4.14: CAMx spatial O3 output in the Terrace - Kitimat valley airshed for the Test case on August 15th, 2010 at: a) 08:00, b)10:00, c) 12:00, d) 14:00, e) 16:00, f) 18:00, g) 20:00 and h) 22:00 PST. The green line outlines the Terrace - Kitimat valley arished.￿, ￿ and ￿ identify the communities Terrace, Kitimat and Hartley Bay respectively, while and identify hotspot locationsLakelse Lake and Miskatla Inlet. Parallel lines identify the locations of roads.93From these figures one can see how on the morning of August 14th a northerly wind (outflow)advected emissions from land over water where photochemistry began to occur around noon(see figure 4.13 (c)). Between 14:00 and 16:00 PST the flow reversed and air began to advectnorthward though elevated mixing ratios were limited to the southern portion of the valley.On August 15th the onshore wind commenced at 12:00 (see figure 4.14 (c)), and the airmassonce again passed over the main emission sources in Kitimat and was adverted further up theTKV where it encounterd elevated biogenic VOC emissions. Given the shallow PBLH therewas a limited volume in which to disperse these pollutants and the result was elevated O3near Lakelse Lake at 16:00 PST (figure 4.14 (e)). The maximum mixing ratio at this locationreached 67.6 ppb, which was a 55.8% increase compared with the control-case emissions.Multiday recirculation has been hypothesized in other areas (Vancouver, Los Angeles, Spainand Israel) and these photochemically-aged air masses are more likely to be sensitive to NOXemissions (Sillman, 2002). Given the nature of the proposed emissions, it is thus not surprisingthat this combination of meteorology and emissions are the ingredients for elevated O3 in theTKVA. While at no time did these elevated mixing ratios exceed the provincial one-hourambient O3 objective, some sources of uncertainty exist from the model results which mayhave led to under-predicting O3 mixing ratios:• In reality temperatures from August 13th through August 16th were much warmer thanWRF-modelled temperatures, exceeding 35 °C on August 14th in the TKV (see appendixC), and• in reality, wind direction had an onshore-offshore diurnal variation from August 14ththrough August 17th, however WRF output only shows this exchange from August 14ththrough August 15th.Based on this information bVOC emissions would have been even higher than those output byMEGAN. Also, horizontal recirculation has been shown to lead to increased O3 air pollution(Steyn, 2003). The combination of extra bVOC emissions with additional days of recirculationmay have led to even greater O3 production during this time. Additional research is warranted94to investigate whether this would have happened.Historical Frequency of Meteorological Factors Which Contribute to Elevated O3Two of the the above ingredients were parameters measured at meteorological monitoringstations in Kitimat, namely temperature and wind direction. Monitoring data from January2001 until August 2015 was reviewed and dates meeting the criteria of two or more consecutivedays with daily maximum temperatures greater than 29 °C and onshore-offshore wind patternsare listed in table 4.1. It is not possible to know what the sea surface temperature was onthose days as well as the height of the PBLH, as those were not measured.Table 4.1: Dates meeting screening criteria for elevated O3 in the Terrace - Kitimat valleyYear Period Event length (number ofconsecutive days)20012002 Jun 12 - 14 32003 July 28 - 30 32004 Jun 16 - 23 8July 13 - 15 3Aug 13 - 15 32005 Aug 8 - 13 72006 Jun 9 - 13 5July 21 - 22 22007 July 11 - 12 220082009 June 2 - 6 5July 25 - Aug 6 132010 July 7 - 8 2Aug 12 - 17 620112012 Aug 15 - 17 32013 July 28 - Aug 1 5Sep 13 - 14 22014 July 11 - 14 4Aug 11 - 13 32015 July 4 - 9 6July 17 - 19 3The results from this brief analysis suggest that the median frequency meeting the abovecriteria is 1 occurrence per year with a median length of 3 days. The previous three summershave each met this criteria twice, with a median length of 3.5 days (mean of 3.8 days). From95this information is is plausible to assume that elevated O3 ingredients present themselvesannually. More research is needed to confirm the PBLH for those periods, as well confirmwhether there was a differential in sea and land surface temperature.4.2.4 Effect of a Wildfire on Ozone ProductionChanging the emissions inventory in SMOKE is not the only method to modify emissionsinput to CAMx. The other method, used in this section, is modification of the model’s ICs.Modifying the ICs provides a one-time opportunity to change the model outcomes, potentiallyinfluencing up to 84 hours of the model run.The rationale for the importance of investigating ICs is provided in this section. Only indus-trial, mobile and biogenic emissions were included in the emission inventory though in realitythere were a host of uncounted emissions emitted into the model domains during both peri-ods. The addition of low-level area sources was intended to account for any and all other NOXemissions within Kitimat itself, though other sources were not represented. One of these wasa forest fire within the 4 km model domain which had elevated emissions on August 14th -16th, 2010 and led to high concentrations of PM2.5 in Smithers on those days. Images takenfrom a viewpoint near Smithers on August 14th and 15th (see figures 4.15 (a) and 4.15 (b))illustrate the polluted nature of the BV.The forest fire’s emissions were not captured in the SMOKE emissions inventory, and, as theyoccurred within the domain, they were also not captured as part of the model’s boundaryconditions. The only method to investigate their effect was to initiate the model on the datesleading up to August 15th and analyze the output. Changes to resulting ambient concentrationscould then be attributed to the forest fire.Method for Incorporating Forest Fire Emissions into the CAMx Model RunsAs mentioned in section 2.4, boundary and ICs were provided to CAMx from MOZART.This model relies on a variety of emission inputs including anthropogenic emissions based on96Figure 4.15: Haze in the Bulkley valley from a nearby forest fire on August 14th (left) andAugust 15th (right), 2010. PM2.5 concentrations: 84 µg/m3 and 157 µg/m3 respectively. Photocredit: Eric Piercethe ARCTAS global emissions inventory (http://bio.cgrer.uiowa.edu/arctas/emission.html), fire emissions from FINN-v1 (Wiedinmyer et al., 2011) as well as biogenic emissionsdeveloped using MEGAN. Fire emissions from FINN-v1 include emissions of trace gas andparticles from wildfires, agricultural fires and prescribed burning but exclude biofuel use andgarbage burning (Wiedinmyer et al., 2011). Emissions are satellite derived, based on observa-tions from MODIS instruments onboard the NASA Terra and Aqua polar-orbiting satellitesand have a global resolution of 1 km2. Forest fire information is available for download athttp://bai.acom.ucar.edu/Data/fire/. Fire emissions which occurred within the variousdomains are illustrated in figure 4.16.IC images for NO, NO2, and O3 are presented in figures 4.17, 4.18 and 4.19 respectively forJuly 29th and August 12th - 14th. From these images it is possible to see how forest fireemissions captured by FINN-v1 change ICs in MOZART and ultimately those of CAMx.97Figure 4.16: Forest fire locations within the 12 km model domain on August 14th 2010.indicate forest fire locations. The red line outlines the CAMx 12 km domain, the green lineoutlines the Terrace - Kitimat valley airshed and the purple line outlines the Bulkley valleyairshed. ￿, ￿ and ￿ identify the locations of Terrace, Kitimat and Smithers respectively.98Figure 4.17: Spatial NO fields generated by the Model of Ozone and Related Chemical Tracers for use as initial conditions to CAMxfor: a) July 29th, b) August 12th, c) August 13th and d) August 14th. Forest fire locations are shown on August 14th. The redline outlines the CAMx 12 km domain, the green line outlines the Terrace - Kitimat valley airshed and the purple line outlines theBulkley valley airshed. ￿, ￿ and ￿ identify the locations of Terrace, Kitimat and Smithers respectively.99Figure 4.18: Spatial NO2 fields generated by the Model of Ozone and Related Chemical Tracers for use as initial conditions toCAMx for: a) July 29th, b) August 12th, c) August 13th and d) August 14th. Forest fire locations are shown on August 14th. Thered line outlines the CAMx 12 km domain, the green line outlines the Terrace - Kitimat valley airshed and the purple line outlinesthe Bulkley valley airshed. ￿, ￿ and ￿ identify the locations of Terrace, Kitimat and Smithers respectively.Figure 4.19: Spatial O3 fields generated by the Model of Ozone and Related Chemical Tracers for use as initial conditions to CAMxfor: a) July 29th, b) August 12th, c) August 13th and d) August 14th. Forest fire locations are shown on August 14th. The redline outlines the CAMx 12 km domain, the green line outlines the Terrace - Kitimat valley airshed and the purple line outlines theBulkley valley airshed. ￿, ￿ and ￿ identify the locations of Terrace, Kitimat and Smithers respectively.100Temporal ResultsFigure 4.20 shows how initializing CAMx using ICs from August 12th - 14th changes ambientO3 at Lakelse Lake, the location with the highest mixing ratio only for the model initializedon August 14th. O3 mixing ratios for the model run initialized on August 12th are identicalto those from the July 29th initialization after only 24 hours. It takes approximately 48 hoursfor this to happen to the run initialized on August 13th. Neither of these runs incorporate theforest fire’s emissions, and resulting O3 mixing ratios on August 15th are identical to thosefrom the July 29th initialization. MOZART’s inclusion of forest fire emissions in the ICs onAugust 14th adds a further 10 ppb O3, leading to a peak mixing ratio just under 80 ppb.Differences in ambient mixing ratios persist for an additional 24 hours with final convergenceto the original time series approximately 84 hours after initialization. It should be notedthat the provincial air quality advisory threshold for O3 is a one-hour mixing ratio of 82 ppb.Therefore, the combination of industrial emissions with the ingredients required for elevatedO3 and wildfires came close to, but did not exceed this provincial objective.Figure 4.20: O3 time series at Lakelse Lake with varied initial conditions. Test Case time seriesinitialized on July 29th, August 12th, August 13th and August 14th are shown as dashed red,purple, blue and turquoise lines respectively. Control case time series for all dates are shownas dotted black lines. Initial O3 mixing ratios for each run are shown as coloured triangles.101Spatial ResultsCAMx output suggests that for a case such as a wide-spread forest fire, increased O3 wouldnot be limited to the point of maximum mixing ratio. Indeed, O3 mixing ratios were increasedacross the entire model domain. Figure 4.21 shows ambient O3 at 16:00 PST on August 15thfor the test case initiated on July 29th and also for the Test case initiated on August 14th.102Figure 4.21: CAMx spatial O3 output for the Terrace - Kitimat valley airshed on August 15th, 2010 at 16:00 PST for Test Casemodel initialization on: a) July 29th and b) August 14th. The green line outlines the Terrace - Kitimat valley arished. ￿, ￿ and￿ identify the communities Terrace, Kitimat and Hartley Bay respectively, while and identify hotspot locations Lakelse Lakeand Miskatla Inlet. Parallel lines identify the locations of roads.103Chapter 5Results and Discussion III: AnsweringQuestion 3The cycling of reactions 1.5 through 1.9 terminates through either reaction 1.11 or 1.12.Based on which terminal reaction dominates it is possible to determine the sensitivity ofO3 production to emissions of NOX or VOCs. In low-NOX environments 1.11 dominatesand O3 production varies linearly with increasing NO, while in high-NOX environments 1.12dominates and O3 production varies linearly with increasing VOCs but inversely to increasingNO2 (Jacob, 1999).These sensitivities are applicable both spatially and temporally. For example, in urban centres(or large NOX producing areas) production of O3 is generally VOC-limited while in less urbanareas O3 production is limited by emissions of NOX (Jacob, 1999). Locations which are NOX -sensitive at one time of day can be VOC-sensitive at other times of the day (Sillman and West,2009), and NOX -sensitive locations in the summer can transition to be VOC-sensitive in theautumn (Jacob et al., 1995). The final question in this research is aimed at understanding thepresent and futures sensitivity of O3 in the TKVA to emissions of NOX and VOCs:• What is the current and future O3 sensitivity of the Terrace - Kitimat valley airshed to104emissions of NOX and VOCs?As management strategies for O3 must take into account the complexities and non-linear natureof its production, this question is often asked when policy makers are attempting to maximizeO3 reduction strategies (Sillman and He (2002), Jacob (1999)). In this application NOX -VOCsensitivity is assessed to better understand the future state of the TKV and whether or notfurther increases in NOX would lead to additional O3.5.1 Method for Determining NOX-VOC SensitivityOnly model output from the summer period was used to assess the sensitivity of O3 formationin the TKVA to emissions of NOX and VOCs. The spring period was omitted as O3 was notgenerated in meaningful quantities along the valley bottom. Sensitivity designations are oftenbased on modelled change in O3 relative to decreases in NOX or VOCs from various emissionscenarios (Sillman and He, 2002). An important parameter for this analysis is total reactivenitrogen (NOy), which includes the gases in NOX as well as NO3, HNO3, HONO, N2O5,HO2NO2, peroxyacetyl nitrate (PAN) and organic nitrates, but excludes NH3 (Neuschuler,D., 2006). Ainslie et al. (2013) demonstrated how, in the LFV, the ratio of [O3]/[NOy] ≈ 7can be used to differentiate NOX and VOC-sensitive areas, where [O3]/[NOy] < 7 is sensitiveto VOCs and [O3]/[NOy] > 7 is sensitive to NOX . This result is similar to a result from(Castell et al., 2009) who found that 6.11 was an appropriate ratio considering an O3 responseof at least 3 ppb to changes in NOX emissions in an airshed with a natural gas power plantemitting 4.45 tonnes of NOX per day. Sillman and He (2002) proposed a range of 11-15 as thetransition range for areas with O3 mixing ratios less than 80 ppb. Given the above findings,this research applies a broad transition range of 7 - 15 for the TKVA; one area for furtherstudy could be to reduce this range based on localized findings and further sensitivity tests.In order to facilitate this analysis, a transect of the TKV was created. This transect follows themain trajectory of the large industrial plumes, both to the north (when the wind is onshore),and to the south (when the wind is offshore). The domain was also divided into north and105south sections. This enabled separate analysis of each section, depending on the dominantwind direction. Figure 5.1 shows the transect and north-south divide.Figure 5.1: The north - south transect and north-south divide of the Terrace - Kitimat valley.The orange line indicates the transect and the yellow line identifies the divide. The green lineoutlines the Terrace - Kitimat valley arished. ￿,￿ and￿ identify the communities Terrace,Kitimat and Hartley Bay respectively, while and identify hotspot locations Lakelse Lakeand Miskatla Inlet.For each half of the domain, the [O3]/[NOy] ratio was extracted for the hours where thedifference between Control and Test Case modelled O3 was 4 ppb or greater. This enabledthe analysis of only those hours where O3 was produced. These differences were calculatedat Lakelse Lake (to the north) and Miskatla Inlet (to the south). In other studies the hoursselected have been in the afternoon, typically between 1300 and 1600 local time (Steyn etal. (2012), Stein et al. (2005), Torres-Jardon and Keener (2006)), though in the TKVA O3production often commenced earlier (noon) and continued later (until 1800). Thus it was morerelevant to use hours where the difference was 4 ppb or greater rather than a set time. The106dates and hours used for analysis of both the northern and southern halves of the domain arelisted in table 5.1. As the dominant wind direction was northerly for this period, there wereonly two days that satisfied the criteria for the northern half. (See the time series of controland test cases for Lakelse Lake and Miskatla Inlet in figure 4.5)Table 5.1: Hours used for [O3]/[NOy] analysis.South NorthDate Time Range Number of Hours Date Time Range Number of HoursJuly 31 14:00 - 16:59 3 Aug 15 13:00 - 19:59 7Aug 3 14:00 - 19:59 6 Aug 16 10:00 - 14:59 5Aug 4 10:00 - 16:59 7Aug 12 13:00 - 16:59 4Aug 13 13:00 - 19:59 7Aug 14 11:00 - 21:59 11Aug 15 08:00 - 15:59 85.2 Sensitivity ResultsThe mean and median [O3]/[NOy] ratios were calculated for each half of the domain, along withthe 5th and 95th percentiles. Sillman and He (2002) identify the transition zone between VOCand NOX -sensitive regions to be the area between the 95th percentile of the VOC-sensitivearea and the 5th percentile of the NOX sensitive area. In this analysis the 7 - 15 range isused because of its breadth. Transect results from both halves were stitched together and arepresented in figure 5.2 for both the Control and Test cases. Also presented are box-whiskerplots at five locations along the TKV: Hartley Bay, Miskatla Inlet, Kitimat, Lakelse Lake andTerrace.As can be seen from these figures, the full construction of all industrial projects in the TKVAwould result in a change in sensitivity of large sections of the TKV on O3-producing daysfrom NOX to VOC-sensitive, particularly from Miskatla Inlet to 10 km south of Lakelse Lake.The transect passing through Kitimat shows the median [O3]/[NOy] ratio decreasing from11 to 1, indicating that without additional VOC emissions, model output represents a worstcase representation of O3 concentrations in this area. At Lakelse Lake the median [O3]/[NOy]107Figure 5.2: Median, mean, 5th and 95th percentile [O3]/[NOy] ratios along the Terrace -Kitimat valley transect defined in figure 5.1 during hours identified in table 5.1 for the: a)Control and b) Test cases. Red line is the median while black line is the mean. Light shadingindicates ares where [O3]/[NOy] < 7 (sensitive to emissions of VOCs), violet shading indicatesareas where 7 < [O3]/[NOy] < 15, and purple shading indicates areas where [O3]/[NOy] > 15(sensitive to emissions of NOX). The vertical yellow line represents the north-south divide.Box plots showing the interquartile range and outliers for transect cells matching the samelatitude as Hartley Bay, Miskatla Inlet, Kitimat, Lakelse Lake and Terrace are shown in green,cyan, blue purple and yellow shading respectively.ratio decreases from 45 to 12, suggesting that in this area additional NOX would still lead toadditional O3. Similar to Lakelse Lake, south of Miskatla Inlet to Hartley Bay the sensitivity108also changes from NOX -sensitive to a transitional sensitivity. Note that photochemically agedplumes are increasingly more NOX -sensitive (Sillman, 2002); additional research is neededto investigate the balance between predicted NOX -insensitive (and transitional) areas withNOX -sensitive plumes.Figure 5.3, presents domain-wide renderings of the median [O3]/[NOy] ratio for (a) Controland (b) Test Cases. This figure illustrates how the change in sensitivity from NOX to VOCsis greatest in and around Kitimat (approximately 25 km to the north and 25 km to the south)and while the changes are greatest along the trajectory of the emissions plume, the decreaseis spread over the entire width of the TKV where it is narrow. Areas with [O3]/[NOy] < 7 areinsensitive to NOX emissions, suggesting that further NOX would not lead to increased O3 inthose areas. Purple areas are still NOX -sensitive, indicating that additional NOX would stilllead to increased O3. The violet area is sensitive to both NOX and VOCs though it is unclearhow much O3-generating potential remains in this zone. Additional research is warranted tofurther our understanding of the 7 < [O3]/[NOy] < 15 area, as this zone occupies a largepercentage of the TKV bottom.109Figure 5.3: Median [O3]/[NOy] ratios in the Terrace - Kitimat valley airshed during hoursidentified in table 5.1 for the: a) Control and b) Test cases. Light shading indicates areswhere [O3]/[NOy] < 7 (sensitive to emissions of VOCs), violet shading indicates areas where 7< [O3]/[NOy] < 15 (zone of transition), and purple shading indicates areas where [O3]/[NOy] >15 (sensitive to emissions of NOX). The green line outlines the Terrace - Kitimat valley arishedand ￿, ￿ and ￿ identify the communities Terrace, Kitimat and Hartley Bay respectively,while and identify hotspot locations Lakelse Lake and Miskatla Inlet.110Chapter 6Summary and ConclusionThe WRF, SMOKE, MEGAN, MOZART and CAMx models were run for two periods in2010 in order to assess the potential change in ambient mixing ratios of O3 and its precursorsthat may arise from the construction of large industrial facilities in the TKVA. Control andtest cases were developed and applied for each period, the former for model evaluation andthe latter to assess pollutant change. The current and future O3 sensitivity in the TKVA toemissions of NOX and VOCs was also assessed.6.1 Findings: Research Questions RevisitedThe findings of this research are put in the context of the original research questions:1. Can spring and summer photochemical O3 be replicated by a model in a manner thatis fit for the purpose of investigating worst-case concentrations that may result fromproposed industrial emissions in a constrained coastal airshed with complex terrain,such as the Terrace - Kitimat valley?Photochemical O3 was replicated by a model in a manner that was fit for the purposeof investigating worst-case concentrations albeit not without some difficulties initially.111While not explicitly mentioned in the question, NO and NO2 were included in the anal-ysis because of the important relationship between these gases and O3 as described byreactions 1.1 and 1.2. Model fitness was evaluated based on a combination of the overallmodelled and observed time series, scatterplots and five statistical metrics (NMB, MB,RMSE, MAE and r). Daytime O3 maxima were replicated reasonably however the modelwas originally unable to recreate overnight O3 titration by NO. NO and NO2 time seriesand statistical metrics were originally poor when compared with observed mixing ratios.Because of this the original CAMx Control Case output was deemed unfit. Sensitivitytests were developed to investigate the cause of the lack of overnight O3 titration. It wasdeemed that the original emissions inventory was inadequate to characterize overnightO3, NO and NO2 mixing ratios. Results improved with the addition of low-level NOXsources to account for sources absent in the emissions inventory, as these sources allowedovernight titration to occur without affecting daytime O3 production. The model wasthen determined to be fit for the purpose of investigating worst-case O3 concentrationsand thus the approach described in chapter 2 was deemed to be an appropriate modellingapproach.2. Should all proposed industrial facilities in the Terrace - Kitimat valley airshed be con-structed,(a) Where are the locations of the O3 maxima and what are the worst-case concentra-tions that could occur there?(b) What meteorological factors contribute to enhancing or reducing O3 concentrationsin the Terrace - Kitimat valley airshed?Results for the spring period illustrated that, despite the naturally elevated O3 in thespring, the addition of NOX from industrial sources led to O3 titration overnight and alsoduring some daytime hours along low elevations of the TKV. This was most prevalent inthe near-field around Kitimat but occurred upwards of 80 km downwind, occasionally inTerrace. The meteorological conditions which led to the greatest change (i.e.: decrease)112in springtime O3 were low PBLHs and consecutive days of recirculating wind. ModestO3 production was found during afternoon hours outside of this main plume and onvalley walls.Results from the summer period illustrated how the addition of O3 precursors contributedto an overall increase in O3 downwind of Kitimat. O3 production varied with winddirection. To the south the hotspot was determined to be Miskatla Inlet, which actedas a trap for pollutants returning northward after advecting southward earlier in themorning hours. To the north the hotspot was Lakelse Lake, a location surrounded byvegetation that emitted high quantities of isoprene when temperatures were greater than29 °C. These findings are in general agreement with those of Castell et al. (2010), whofound O3 maxima between 40 and 140 km away from the hypothetical emission source,and who also found that increases were greatest in NOX -sensitive areas, i.e.: locationswith high bVOC and low NOX emissions.In the summer the meteorological conditions which led to the greatest change (i.e.: in-crease) in O3 were high temperatures (these led to high bVOC emissions), low PBLHs,differential heating of the land and ocean surface temperatures and consecutive days ofrecirculating wind. These conditions led to a greater than 55% increase in O3 produc-tion at both Miskatla Inlet and Lakesle Lake. The frequency of these meteorologicalconditions was estimated to be one per year with a median length of three days based onan analysis of meteorological monitoring data in Kitimat over the past 15 years. Whenemissions from a nearby forest fire were included in the emissions inventory (through theICs input file) the maximum O3 at these locations increased by 10 ppb and approachedthe provincial threshold for O3-based air quality advisories.3. What is the current and future O3 sensitivity of the Terrace - Kitimat valley airshed toemissions of NOX and VOCs?The modelled [O3]/[NOy] ratio during hours conducive to photochemistry was used todetermine the O3 sensitivity of the TKVA. The TKVA is currently very sensitive to NOX113emissions however the full construction of all proposed industrial projects would likelychange the sensitivity of the large portions of the valley. From Miskatla Inlet to 10 kmsouth of Lakelse Lake the change would be from NOX to VOC-sensitive (especially inand around Kitimat). South of Miskatla Inlet and in the area surrounding Lakelese Lakethe change would be from NOX -sensitive to a transitional sensitivity. In VOC-sensitiveareas further NOX emissions would not likely lead to increased O3, though in the thezone of transition there may still be some O3 production potential.6.2 Future StudyAn interesting outcome of this research was the determination of the ingredients of elevatedO3 production in the TKV. It is startling that the ingredients are almost the same for springand summer periods, yet the predicted outcomes are vastly different, based on higher sum-mer temperatures leading to increased biogenic emissions. Some confirmation of this seemsrequired, i.e. can high biogenic emissions in the TKV be confirmed (perhaps through anothermodel, an evaluation of landcover, or (even better) through monitoring). An ideal location forthis would be Lakelse Lake where there is currently a wet deposition station operated by thealuminum smelter. This seems prudent as measured temperatures during the summer modelperiod exceeded those from WRF output and it is possible that bVOCs emissions were evenhigher.Reducing uncertainty related to model inputs is also warranted. With respect to meteorologicalinputs, the RF statistical model identified RH to be an important predictor for O3 in the TKV,yet it was also the parameter that performed poorest in the WRF evaluation. It is unclearhow this affected model results though improvements to modelled humidity seem warranted.Also, improving WRF temperature output should be attempted, potentially using the nudgingtechnique in WRF’s initialization.With respect to the emissions inputs, a thorough emissions inventory is warranted for futuremodelling studies. As many local emission sources as possible should be included, even though114the time required to complete this task could be daunting. NPRI data should be supplementedwith other local data where possible.Future modelling studies would benefit by taking a focused approach and only investigatingthose times identified as meeting the criteria for elevated O3. Modelling the spring perioddemonstrated how the spring O3 peak was insensitive to NOX emissions in this area as VOCemissions are low. The summer provided much more interesting results and future work canfocus on those periods where production is likely, especially extended warm periods withhorizontal recirculation.It must be noted that as of spring 2015, O3 monitoring commenced in Terrace, operated bythe BC MOE. Therefore, future model evaluation could occur within the TKVA and not inneighbouring airsheds.Lastly, scientific study has moved from studying single pollutants to studying multiple pollu-tants. Without an understanding of the total change in atmospheric pollutant concentrations,it is not possible to make final conclusions about the total human and environmental risks fromany combination of these potential new sources. Future iterations may desire to calculate theAir Quality Health Index for different areas in the TKV to begin a true effects assessment, asthat was not the role of this research.115ReferencesAinslie B. Steyn, D., Reuten, C., Jackson, P. (2013). A Retrospective Analysis of Ozone For-mation in the Lower Fraser Valley, British Columbia, Canada. Part II: Influence of EmissionsReductions on Ozone Formation. Atmosphere-Ocean, 51 (2), 170–186.Armstrong, J. (1990). Vancouver Geology. Tech. rep. Visited on 2015-09-29. Geological Associa-tion of Canada. Available from: http://www.gac-cs.ca/publications/VancouverGeology.pdf.Banner, A., MacKenzie, W., Haeussler, S., Thomson, S., Pojar, J., Trowbridge, R. (1993).A field guide to site identification and interpretation for the Prince Rupert Forest Region:Parts 1 and 2. Handbook 26. Victoria: British Columbia Ministry of Forests.Bott, A. (1989). A Positive Definite Advection Scheme Obtained by Nonlinear Renormalizationof the Advective Fluxes. Monthly Weather Review, 117, 1006–1015.British Columbia Ministry of Environment (2014). Provincial Air Quality Objective Informa-tion Sheet. Tech. rep. Visited on 2015-09-29. Available from: http://www.bcairquality.ca/reports/pdfs/aqotable.pdf.British Columbia Ministry of Forests (1998a). The ecology of the alpine tundra zone. Tech. rep.Visited on 2015-09-29. Available from: https://www.for.gov.bc.ca/hfd/pubs/docs/bro/bro56.pdf.British Columbia Ministry of Forests (1998b). The ecology of the engelmann spruce - subalpinefir zone. Tech. rep. Visited on 2015-09-29. Available from: https://www.for.gov.bc.ca/hfd/pubs/docs/bro/bro55.pdf.116British Columbia Ministry of Forests (1998c). The ecology of the sub-boreal spruce zone. Tech.rep. Visited on 2015-09-29. Available from: https://www.for.gov.bc.ca/hfd/pubs/docs/bro/bro53.pdf.British Columbia Ministry of Natural Gas Development (2015a). Everything you want to knowabout LNG in BC. Visited on 2015-09-29. Available from: https://engage.gov.bc.ca/lnginbc/b-c-s-lng-story/#where-is-it-coming-from.British Columbia Ministry of Natural Gas Development (2015b). Explore BC’s LNG projects.Visited on 2015-09-29. Available from: https : / / engage . gov . bc . ca / lnginbc / lng -projects/.Canadian Council of Ministers of the Environment (2012a). Guidance Document on Achieve-ment Determination, Canadian Ambient Air Quality Standards for Particulate Matter andOzone. Tech. rep. Visited on 2015-09-29. Available from: http://www.ccme.ca/files/Resources/air/aqms/pn_1483_gdad_eng.pdf.Canadian Council of Ministers of the Environment (2012b). Guidance Document on Air ZoneManagement. Tech. rep. Visited on 2015-09-29. Available from: http://www.ccme.ca/files/Resources/air/aqms/pn_1481_gdazm_e.pdf.Canadian Council of Ministers of the Environment (2012c). The Air Quality ManagementSystem: Federal, Provincial and Territorial Roles and Responsibilities. Tech. rep. Visited on2015-09-29. Available from: http://www.ccme.ca/files/Resources/air/aqms/pn_1475_roles_and_respn_final_e.pdf.Canadian Press (2013). Greenhouse gas emissions from LNG defended by B.C. premier. Visitedon 2015-09-29. Available from: http://www.cbc.ca/news/canada/british-columbia/greenhouse-gas-emissions-from-lng-defended-by-b-c-premier-1.2424982.117Castell, N., Stein, A., Mantilla, E., Salvador, R., Millan, M. (2009). Evaluation of the use ofphotochemical indicators to assess ozone-NOX -VOC sensitivity in the Southwestern IberianPeninsula. J Atmospheric Chemistry, 63, 73–91.Castell, N., Mantilla, E., Salvador, R., Stein, A., Millan, M. (2010). Photochemical modelevaluation of the surface ozone impact of a power plant in a heavily industrialized area ofsouthwestern Spain. J Environ. Manage. 91 (3), 662–676.Chang, J., Hanna, S. (2004). Air quality model performance evaluation. J Meteorol. Atmos.Phys. 87, 167–196.Chin, M., Jacob, D., Munger, J., Parrish, D., Doddridge, B. (1994). Relationship of ozone andcarbon monoxide over North America. J Geophysical Research. 99 (D7), 14565–14573.Colella, P., Woodward, P. (1984). The Piecewise Parabolic Method (PPM) for Gas- dynamicalSimulations. J Computational Physics, 54, 174–201.Crutzen, P. (1979). The role of NO and NO2 in the chemistry of the troposphere and strato-sphere. Ann. Rev. Earth Planet Sci. 7, 443–472.Crutzen, P. (1995). Overview of Tropospheric Chemistry: Developments during the Past Quar-ter Century and a Look Ahead. Faraday Discuss. 100, 1–21.Demarchi, D. (2010). Coast Mountains Ecoprovince. Visited on 2015-09-29. Available from:http://tinyurl.com/au3kfu.Dennis, R., Fox, T., Fuentes, M., Gilliland, A., Hanna, S., Hogrefe, C., Irwin, J., Rao, S.,Scheffe, R., Schere, K., Steyn, D., Venkatram, A. (2010). A framework for evaluating regional-scale numerical photochemical modelling systems. Environ. Fluid Mech. 10, 471–489.Derwent, R. (1999). Reactive Hydrocarbons and Photochemical Air Pollution. Reactive Hy-drocarbons in the Atmosphere. Ed. by C. Hewitt. London: Academic Press.118Douglas, T. (2012). Integrating ecosystem restoration into forest management: practical exam-ples for foresters. Tech. rep. Visited on 2015-09-29. Available from: http://chapter.ser.org/westerncanada/files/2012/08/2003-Integrating-Ecosystem-Restoration-into-Forest-Management.pdf.Emmons, L., Walters, S., Hess, P., Lamarque, J.-F., Pfister, G., Fillmore, D., Granier, C.,Guenther, A., Kinnison, T., Laepple, T., Orlando, J., Tie, X., Tyndall, G., Wiedinmyer, C.,Baughcum, S., Kloster, S. (2010). Description and evaluation of the Model for Ozone andRelated chemical Tracers, version 4 (MOZART-4). Geosci. Model Dev. 3, 43–67.ENVIRON (2013). User’s Guide: Comprehensive Air Quality Model with Extensions, CAMx(Version 6.0). ENVIRON International Corporation. Novato, CA, USA.ESSA-Technologies, Laurence, J., Limnotek, International, R. S., Alcan, R. T., University,T., Consultants, T., Illinois, U. of (2013). Sulphur Dioxide Technical Assessment Report inSupport of the 2013 Application to Amend the P2-00001 Multimedia Permit. Final TechnicalReport Volume 2. Presented to the BC Ministry of Environment.ESSA-Technologies, Laurence, J., International, R. S., University, T., Consultants, T., Illi-nois, U. of (2014). Kitimat Airshed Effects Assessment. Final Report. Prepared for the BCMinistry of Environment.ESWG (1995). A National Framework for Canada. Tech. rep. Agriculture and Agri-FoodCanada, Research Branch, Centre for Land and Biological Resources Research and Environ-ment Canada, State of Environmental Directorate, Ecozone Analysis Branch. Ottawa/Hull:Ecological Stratification Working Group.Fall, R. (1999). Biogenic Emissions of Volatile Organic Compounds from Higher Plants. Re-active Hydrocarbons in the Atmosphere. Ed. by C. Hewitt. London: Academic Press.Fenger, J. (2002). Urban air quality. Air Pollution Science for the 21st Century. Ed. by J.Austin, P. Brimblecombe, W. Sturges. Oxford: Elsevier.119Finlayson-Pitts, B., Pitts Jr, J. (2000). Chemistry of the Upper and Lower Troposphere.Theory, Experiments and Applications. San Diego, CA, USA: Academic Press.Groupe International Des Importateurs de Gaz Naturel Liquifie (2009). Information PaperNo. 1 - Basic Properties of LNG. Tech. rep. Visited on 2015-09-29. Available from: http://www.kosancrisplant.com/media/5648/1-lng_basics_82809_final_hq.pdf.Guenther, A., Zimmerman, P., Harley, P. (1993). Isoprene and Monoterpene Emission RateVariability: Model Evaluations and Sensitivity Analyses. J Geophysical Research, 98 (D7),12609–12617.Hastie, T., Tibshirani, R., Friedman, J. (2009). The Elements of Statistical Learning. DataMining, Inference and Prediction. New York, NY, USA: Springer.Henolson, A. (2014). Personal correspondence with Trinity Consultants.Hoffman, E. (2015). Personal correspondence with the Ministry of Environment.Hsieh, W. (2009). Machine Learning Methods in the Environmental Sciences. Neural Networksand Kernels. Cambridge, UK: Cambridge University Press.Jacob, D. (1999). An Introduction to Atmospheric Chemistry. Princeton, NJ, USA: PrincetonUniversity Press.Jacob, D., Horowitz, L., Munger, W., Heikes, B., Dickerson, R., Artz, R., Keene, W. (1995).Seasonal transition from NOx - to hydrocarbon-limited conditions for ozone production overthe eastern United States in September . J Geophysical Research, 100 (D5), 9315–9324.Jenkin, M., Clemitshavm, K. (2002). Ozone and other secondary photochemical pollutants.Air Pollution Science for the 21st Century. Ed. by J. Austin, P. Brimblecombe, W. Sturges.Oxford: Elsevier.120Karamchandani, P., Vijayaraghavan, K., Yarwood, G. (2011). Sub-Grid Scale Plume Mod-elling. Atmosphere, 2, 389–406.Krupa, S. (1997). Air Pollution, People and Plants. St. Paul, MN, USA: APS Press.Meidinger, D., Pojar, J. (1991). Ecosystems of British Columbia. Tech. rep. Special ReportSeries 6. Victoria: British Columbia Ministry of Forests.Monks, P. S. (2000). A review of the observations and origins of the spring ozone maximum.Atmospheric Environment, 34, 3545–3561.National Oceanic & Atmospheric Administration (2014). Weather Research & ForecastingWRF Model. Visited on 2015-09-29. Available from: http://www.atdd.noaa.gov/?q=node/41.Neuschuler, D. (2006). Air Monitoring Instrumentation: Nitrogen Oxides (NOy). Visited on2015-09-29. Available from: http://www3.epa.gov/ttnamti1/files/2006conference/neuschuler.pdf.Otte, T., Pleim, J. (2010). The Meteorology-Chemistry Interface Processor (MCIP) for theCMAQ modeling system: updates through MCIPv3.4.1. Geosci. Model Dev. 3, 243–256.Ozone Secretariat (2015). Treaties and Decisions. Visited on 2015-09-29. Available from: http://ozone.unep.org/en/treaties-and-decisions.Pierce, E. (2015). Personal correspondence with the Ministry of Environment.Reid, P., Schroeder, M., Spagnol, J., Gallagher, J., Wang, J. (2014). Air Quality TechnicalData Report: LNG Canada Export Terminal. Tech. rep. Calgary: Stantec Consulting Ltd.Sharkey, T., Wiberley, A., Donohue, A. (2008). Isoprene Emission from Plants: Why and How.Annals of Botany, 101, 5–18.121Sillman, S. (2002). The relation between ozone, NOX and hydrocarbons in urban and pol-luted rural environments. Air Pollution Science for the 21st Century. Ed. by J. Austin, P.Brimblecombe, W. Sturges. Oxford: Elsevier.Sillman, S. (2012). Overview: Tropospheric Ozone, Smog And Ozone-NOX-VOC Sensitivity.Tech. rep. Visited on 2015-09-29. Available from: http://www- personal.umich.edu/~sillman/Sillman-webOZONE.pdf.Sillman, S., He, D. (2002). Some theoretical results concerning O3-NOX -VOC chemistry andNOX -VOC indicators. J Geophysical Research, 107 (D22), 26–1–26–15.Sillman, S., West, J. (2009). Reactive nitrogen in Mexico City and its relation to ozone-precursor sensitivity: results from photochemical modls. J Atmos. Chem. Phys. 9, 3477–3489.Simon, H., Baker, K., Phillips, S. (2012). Complication and interpretation of photochemicalmodel performance statistics published between 2006 and 2012. Atmospheric Environment,61, 124–139.Song, J., Lei, W., Bei, N., Zacala, M., Foy, B. de, Volkamer, R., Cardenas, B., Zheng, J., Zhang,R., Molina, L. (2010). Ozone response to emission changes: a modelling study during theMCMA-2006/MILAGRO campaign. J Atmos. Chem. Phys. 10, 3827–3846.Statistics Canada (2014a). Census for Kitimat, British Columbia. Visited on 2015-09-29. Avail-able from: https://www12.statcan.gc.ca/census-recensement/2011/as-sa/fogs-spg/Facts-csd-eng.cfm?LANG=Eng&GK=CSD&GC=5949005.Statistics Canada (2014b). Census for Smithers, British Columbia. Visited on 2015-09-29.Available from: https://www12.statcan.gc.ca/census-recensement/2011/as-sa/fogs-spg/Facts-csd-eng.cfm?LANG=Eng&GK=CSD&GC=5951043.122Statistics Canada (2015a). Census for Terrace, British Columbia. Visited on 2015-09-29. Avail-able from: https://www12.statcan.gc.ca/census-recensement/2011/as-sa/fogs-spg/Facts-cma-eng.cfm?LANG=Eng&GK=CMA&GC=965.Statistics Canada (2015b). Population of census metropolitan areas. Visited on 2015-10-30.Available from: http://www.statcan.gc.ca/tables-tableaux/sum-som/l01/cst01/demo05a-eng.htm.Stein, A., Mantilla, E., Milan, M. (2005). Using measured and modeled indicators to assessozone-NOx-VOC sensitivity in a western Mediterranean coastalenvironment. AtmosphericEnvironment, 39, 7167–7180.Steyn, D., Ainslie, B., Reuten, C., Jackson, P. (2013). A Retrospective Analysis of OzoneFormation in the Lower Fraser Valley, British Columbia, Canada. Part I: Dynamical ModelEvaluation. Atmosphere-Ocean, 51 (2), 153–169.Steyn, D. (2003). Mesoscale circulations and regional air pollution. Air pollution processes inregional scale. Ed. by D. Melas, D. Syrakov. Springer.Steyn, D., Builtjes, P., Schaap, M., Yarwood, G. (2012). Regional Air Quality Modelling: NorthAmerican and European Perspectives. Environmental Manager, July, 6–8.Stull, R. (2014). Practical Meteorology: An Algebra-based Survey of Atmospheric Science.Available online at http://www.eos.ubc.ca/books/Practical_Meteorology/.Swall, J., Foley, K. (2009). The impact of spatial correlation and incommensurability on modelevaluation. Atmospheric Environment, 43, 1204–1217.Swedish Environmental Protection Agency (2009). Swedish Pollutant Release and TransferRegister: Nitrogen oxides (NOX). Visited on 2015-09-29. Available from: http://utslappisiffror.naturvardsverket.se/en/Substances/Other-gases/Nitrogen-oxides//.123Tiwary, A., Colls, J. (2010). Air Pollution. Measurement, modelling and mitigation, thirdedition. New York, NY, USA: Routledge.Torres-Jardon, R., Keener, T. (2006). Evaluation of Ozone-Nitrogen Oxides-Volatile OrganicCompound Sensitivity of Cincinnati, Ohio. J Air and Waste Management Association, 56,322–333.UNC (2013). User’s Manual: Sparse Matrix Operator Kernel Emissions, SMOKE (Version3.5.1). University of North Carolina at Chapel Hill. Chapel Hill, NC, USA.U.S. Environmental Protection Agency (2015a). Overview of Greenhouse Gases. Visited on2015-09-29. Available from: http://epa.gov/climatechange/ghgemissions/gases/ch4.html.U.S. Environmental Protection Agency (2015b). Photochemical Modeling. Tech. rep. Visited on2015-09-29. Available from: https://www.for.gov.bc.ca/hfd/pubs/docs/bro/bro31.pdf.Vallero, D. (2014). Fundamentals of Air Pollution. Fifth ed. San Diego, CA, USA: AcademicPress.Vingarzan, R. (2013). Personal correspondence with Environment Canada.Wiedinmyer, C., Akagi, S., Yokelson, R., Emmons, L., Al-Saadi, J., Orlando, J., Soja, A.(2011). The Fire INventory from NCAR (FINN): a high resolution global model to estimatethe emissions from open burning. Geosci. Model Dev. 4, 625–641.Willmott, C. (1982). Some comments on the evaluation of model performance. Bulletin of theAmerican Meteorological Society, 63, 1309–1313.Willmott, C., Matsuura, K. (2005). Advantages of the mean absolute error (MAE) over theroot mean square error (RMSE) in assessing average model performance. Climate Research,30, 79–82.124Wilson, G. (2014). Personal correspondence with CAMx Development Team, ENVIRON In-ternational Corporation.WSU (2012). User’s Guide: Model of Emissions of Gases and Aerosols from Nature, MEGAN(Version 2.10). Washington State University. Pullman, WA, USA.Yarwood, G., Yocke, M., Whitten, G. (2005). Updates to the Carbon Bond chemical mechanism:CB05. Final Report 415.899.0703. Novato, CA: Prepared for U.S. E.P.A.Yarwood, G., Jung, J., Whitten, G., Heo, G., Melberg, J., Estes, M. (2010). Updates to the Car-bon Bond mechanism for version 6 (CB6). Presented at the 9th Annual CMAS Conference.Novato, CA 94998: ENVIRON International Corporation.Zhang, L., Brook, J., Vet, R. (2003). A revised parameterization for gaseous dry deposition inair-quality models. J Atmos. Chem. Phys. 3, 2067–2082.125Appendix ARandom Forest Regression Method forModelling Ozone in KitimatThis chapter contains a brief description of the regression model that was used to select datesfor the spring and summer model periods.A.1 BackgroundTwo periods were selected for this research, one in the spring and one in the summer. Aspring period was chosen because O3 is naturally elevated during this time (Monks, 2000) andthe result of adding additional NOX into the TKV warranted some investigation. A summerperiod was selected as summer is a more traditional time of year for O3 episodes due toincreased production of the OH radical through reactions 1.3 and 1.4 (Jacob, 1999). Becausecomputational resources are expensive, running CAMx for the entire spring and summer wasnot possible. Instead, finding a range of dates was desired where either:(a) O3 was already elevated, or(b) conditions conducive for photochemical O3 production were present.126A.2 MethodsA statistical regression model was used to determine, in a scoping capacity, date ranges thatmet the above criteria. As O3 data is limited to two short periods in Kitimat (Sept 16th -Dec 1st, 2010 and May 20th - Nov. 21st, 2011), the model was developed (trained and tested)using air quality and meteorological observations from Smithers. Input predictors originallyincluded: temperature, RH, wind speed, NO, NO2, CO and PM2.5.Three different regression models were tested for both spring and summer to investigate whichcould best emulate an O3 time series using the above predictors.1. multiple linear regression (MLR),2. neural networks (NN), and3. RF (ensemble of Classification and Regression Trees).Background information on all these models can be found in Hsieh (2009) and Hastie et al.(2009).As complete time series are required for both predictor and response variables, only data from2007, 2010 and 2011 were used for the spring while and 2007, 2009, 2010 and 2011 were usedfor the summer (other years have large data gaps for one or more variable). Predictors wereinitially trained to model O3 for the same year, and regression coefficients were developed.(Note that NN and RF models do not develop regression coefficients; for extreme simplicitythe term is used in the context of NN and RF to describe what those models use to predictresponse variables.) Thus each year constituted its own model. Metrics chosen to evaluatemodel performance and fitness were: RMSE, MAE and r. Equations for these metrics arepresented in appendix B.127A.3 ResultsTables A.1, A.2 and A.3 list the evaluation result statistics for each model type and year basedon its ability to recreate observed O3 for the same year. As can be seen, for all metrics inall years the RF model outperformed the other two model types and it was chosen to be thestatistical model for predicting O3 in Kitimat.Table A.1: Model evaluation statistics for spring and summer periods at Smithers: multiplelinear regression modelSpring SummerModel Year n RMSE MAE r n RMSE MAE r(ppb) (ppb) (ppb) (ppb)2007 2080 8.324 6.722 0.823 1640 4.759 3.589 0.8832009 1570 4.287 3.289 0.9082010 2310 6.969 5.541 0.871 1640 4.196 3.181 0.9152011 2360 6.675 5.179 0.883 1250 4.025 3.138 0.888Table A.2: Model evaluation statistics for spring and summer periods at Smithers: neuralnetwork regression modelSpring SummerModel Year n RMSE MAE r n RMSE MAE r(ppb) (ppb) (ppb) (ppb)2007 2080 8.746 6.831 0.804 1640 4.680 3.288 0.8872009 1570 4.592 3.483 0.8942010 2310 6.720 5.199 0.881 1640 4.507 3.389 0.9022011 2360 6.459 4.977 0.891 1250 4.011 3.035 0.889Table A.3: Model evaluation statistics for spring and summer periods at Smithers: randomforest regression modelSpring SummerModel Year n RMSE MAE r n RMSE MAE r(ppb) (ppb) (ppb) (ppb)2007 2080 6.571 4.884 0.895 1640 3.998 2.796 0.9192009 1570 3.867 2.916 0.9262010 2310 5.832 4.372 0.912 1640 3.916 2.917 0.9272011 2360 5.401 3.969 0.925 1250 3.585 2.687 0.913Using regression coefficients determined by particular models for particular years, the RFmodels were then used to predict O3 from different years as well as O3 in Kitimat for the briefperiod with available monitoring data. An ensemble model was also used. Because not all thepredictor variables available in Smithers were available in Kitimat, the number of predictorswas reduced to only those monitored in Kitimat. These included: temperature, RH, wind128speed and PM2.5. It was determined that while all predictors added value to modelling O3in the spring, the models performed better in the summer using only RH and wind speed asinputs. Of all predictors (in both spring and summer), RH was the most important. This wasdetermined by removing each predictor from the model and measuring the increase in RMSEand MAE, and decrease in r. Model results from the RF prediction of O3 in Kitimat are listedin table A.4.Table A.4: Model evaluation statistics for spring and summer periods at Kitimat using modelcoefficients developed from Smithers: random forest regression modelSpring SummerModel Year RMSE MAE r RMSE MAE r(ppb) (ppb) (ppb) (ppb)2007 6.239 5.067 0.778 5.617 4.577 0.6812009 4.897 3.863 0.7322010 6.562 5.339 0.792 5.830 4.813 0.6692011 6.501 5.222 0.784 6.094 5.067 0.644Ensemble 6.403 4.879 0.802 5.422 4.436 0.695As can be seen from table A.4, for all metrics the summer period outperforms the springperiod. The ensemble RF model outperforms the individual models for the spring periodthough the 2009 RF model outperforms the other years, including the ensemble. The causeof the high performance for the 2009 summer model is not known. Based on these results theensemble model was selected for predicting springtime O3 for 2010, while the 2009 model wasselected to model 2010 summertime O3. The resulting O3 time series for the year 2010 forboth spring and summer model periods are shown in figure A.1 (a) and A.1 (b) respectively.The models estimated elevated O3 in 2010 during two periods in the spring (April 14th - 18th,May 2nd - 14th), and two periods in the summer (July 29th - Aug 6th and Aug 12th - 18th).These corresponded to days with low RH, as model results show this to be the most importantpredictor for elevated O3 in the TKV and Smithers. The final model period for each seasonwas selected to start at the beginning of the first case of elevated O3 for each season andcontinue to the end of the second, and are listed in table A.5.129Figure A.1: Modelled O3 time series at Kitimat using RF model coefficients developed inSmithers for: a) spring and b) summer model periods. Final spring model period is from April14th - May 14th, and summer period is from July 29th to August 18th.Table A.5: Dates for the spring and summer model periods (2010).Season Start Date End Date Number of DaysSpring April 14 May 14 31Summer July 29 August 18 21130Appendix BStatistical Metrics Used for Model Eval-uationNote that in each equation Mi represents the ith model output and Oi represents the corre-sponding ith observed output. The difference between Mi and Oi is the error, Ei.Normalized mean biasNMB = 100% ⋅ ∑ (Mi −Oi)∑Oi (B.1)Mean biasMB = 1N⋅￿ (Mi −Oi) (B.2)Root mean squared errorRMSE =￿∑ (Mi −Oi)2N(B.3)Mean absolute errorMAE = 1N⋅￿ ￿Mi −Oi￿ (B.4)131Correlation coefficientr =N∑i=1((Mi − M¯) × (Oi − O¯))￿N∑i=1(Mi − M¯)2 N∑i=1(Oi − O¯)2(B.5)132Appendix CWRF EvaluationC.1 Quantitative EvaluationNote that all figures provided in this section are in Pacific Standard Time (UTC - 8).TemperatureSpring and summer time series of modelled and observed temperature at the three evaluationlocations can be found in figures C.1 (a) and C.1 (b) respectively.These plots show that WRF does a reasonable job modelling temperature. Output is almostalways within the diurnal range measured at each station but is unable to capture the dailypeaks and troughs. Maximum error at all stations occurs around May 9th and August 15th.At Terrace the model temperature is consistently near the warmer end of the observed valueswhereas in Kitimat the modelled temperature is in the middle of the observed values. Fromthese time series it is apparent that model performance is strongest at the airport. Statisticalmetrics are listed in table C.1 below.With the exception of CYXT in the summer period, the overall tendency of WRF is a warm133Figure C.1: Observed and modelled temperature time series at Terrace (top), CYXT (middle)and Kitimat (bottom) for: a) spring and b) summer model periods. Black lines are observedwhile yellow, orange and blue lines are modelled.Table C.1: Modelled - observed temperature statistics from the WRF evaluation.Spring SummerLocation NMB MB RMSE MAE NMB MB RMSE MAE(%) (o C) (°C) (°C) (%) (°C) (°C) (°C)Terrace B.C. Access Centre 6.136 0.601 2.168 1.643 4.317 0.907 2.344 1.788CYXT Terrace Airport 3.980 0.335 2.441 1.896 -0.890 -0.175 2.218 1.687Kitimat Whitesail 4.304 0.362 2.301 1.788 2.538 0.505 2.862 2.209134bias between 0.3 and 0.9 °C. Considering only bias, WRF had the greatest difficulty modellingtemperature in Terrace and performed best at the airport. On the other hand, in the springboth the RMSE and MAE are lowest for Terrace and highest at the airport. Results forKitimat are somewhere in the middle, with the model having a warm bias (0.3 and 0.5 °C forthe spring and summer periods respectively) and in the summer the highest RMSE and MAE.One possible explanation for the model’s inability to capture the highs and lows of the diurnalcycle in Kitimat is that WRF is unable to master the transition from water to land. Indeed,a look at the Landmask variable shows that Kitimat Whitesail is located in a cell adjacent toa water cell (the Kitimat river) even though this feature is 2 km away. This sub-grid scaleeffect is an example of an incommensurability problem, described by Swall and Foley (2009).Relative HumidityRH data is not a direct output from WRF and was calculated by the program NCARCommand Language (NCL) using the temperature, pressure and water vapour mixing ra-tio variables. More information on this process can be obtained from the NCL website:https://www.ncl.ucar.edu/Document/Functions/Built-in/relhum.shtml. (Note that aformula is not given but RH can be determined using guidance from Roland Stull’s onlinetextbook: Practical Meteorology, Chapter 4 : http://www.eos.ubc.ca/books/Practical_Meteorology/Chapters/Ch04-Moist.pdf (Stull, 2014)). In 2010 RH was measured at Terraceand the CYXT airport but not at Kitimat Whitesail. Interestingly, while WRF was able tomodel temperature adequately, as can be seen in figures C.2 (a) and C.2 (b), it had a moredifficult time with RH.The level to which RH was under-modelled in Terrace is quite apparent in both figures of C.2(a). RH troubles are not isolated to the B.C. Access Centre; looking at the airport time series,it is evident that timing of the modelled diurnal RH pattern is not aligned with the timingof the observed RH pattern and appears consistently out of sync by two hours. It is unclearwhat this means from a photochemical perspective, however during the statistical analysis of135Figure C.2: Observed and modelled relative humidity time series at Terrace (top) and CYXT(bottom) for: a) spring and b) summer model periods. Black lines are observed while yellowand orange lines are modelled.O3 production in Smithers high O3 was correlated with low RH values. These issues affect thestatistical metrics, as listed in table C.2.Table C.2: Modelled - observed relative humidity statistics from the WRF evaluation.Spring SummerLocation NMB MB RMSE MAE NMB MB RMSE MAE(%) (%) (%) (%) (%) (%) (%) (%)Terrace B.C. Access Centre -18.165 -11.913 18.912 15.306 -19.168 -11.784 18.144 14.199CYXT Terrace Airport -5.112 -3.458 13.229 10.324 0.891 0.549 12.512 9.978At the B.C. Access Centre in Terrace, RH is under-predicted by approximately 12% in bothspring and summer seasons (MB). This corresponds with an 18.2 and 19.2% NMB for spring136and summer respectively. This is a much poorer showing when compared with the normalizedtemperature biases. The RH RMSE and MAE in Terrace are also quite high, indicating thatthere were some strong discrepancies between modelled and observed values. Despite the timelag, evaluation metrics fare better at the airport. It is interesting that results vary so greatlyover such a short distance given that the airport is less than 5 km from the Access Centre. Itis possible that the difference in elevation between the two stations plays a role in the modeloutput performance (elevation of the airport is over 200 m while the elevation of the AccessCentre is 68 m).PressureIn the TKV pressure is only measured at the CYXT airport. Modelled vs. observed pressuretime series are presented in figures C.3 (a) and C.3 (b). As can be seen, WRF performsoptimally for this parameter at this location.Evaluation metrics for pressure are listed in table C.3. WRF does a satisfactory job modellingatmospheric pressure. Springtime normalized mean bias was less than 0.1% and summertimebias was just over 0.2%. Considering that the these metrics are calculated in Pascals (asopposed to kPa) they stand out as being remarkably low.Table C.3: Modelled - observed pressure statistics from the WRF evaluation.Spring SummerLocation NMB MB RMSE MAE NMB MB RMSE MAE(%) (Pa) (Pa) (Pa) (%) (Pa) (Pa) (Pa)CYXT Terrace Airport 0.081 80.51 191.73 154.42 0.220 218.28 247.31 218.50Wind DirectionFrom an evaluation perspective, wind direction is a very important parameter. Withoutproperly modelling wind direction a model loses its credibility. Usually the U and V windvectors are evaluated as opposed to wind direction as it is not output directly from WRF andmust be calculated from the U and V vectors at 10 m using the following formula:137Figure C.3: Observed and modelled pressure time series at CYXT for: a) spring and b)summer model periods. Black lines are observed and orange lines are modelled.WDIR = atan2(−U, −V ) × ￿180⇡￿ (C.1)Given the north-south orientation of the TKV it was determined that the calculated winddirection could be evaluated as opposed to the U and V vectors themselves. Model andobserved wind direction are presented for the spring and summer periods in figures C.4 (a)and C.4 (b) respectively. From these it is possible to see that WRF does a reasonable job inKitimat and at the CYXT airport but seems to have some difficulty at Terrace. It should benoted that given the dominant north-south wind pattern in this valley, wind directions lessthan 90°were added to 360, thus 405° is equivalent to a northeast wind while 450° is equivalentto an east wind.As is seen in the above figures, WRF’s performance is reasonable, especially in Kitimat where,138Figure C.4: Observed and modelled wind direction time series at Terrace (top), CYXT (mid-dle) and Kitimat (bottom) for: a) spring and b) summer model periods. Black lines areobserved and yellow, orange and blue lines are modelled. Note that to reduce clutter in thefigures, wind from directions 0 - 90° were moved to the 360 - 450° range. Therefore, 405° isequivalent to a northeast wind while 450° is equivalent to an east wind.despite the obvious bias, the diurnal flow is adequately captured. Records of Kitimat White-sail’s wind vane maintenance were recently checked to confirm that observational error wasnot a problem (20° is the approximate difference between true north and magnetic north inKitimat) and in looking at the detailed topography of the location of Kitimat Whitesail it is139possible to explain why, when there is a southerly wind, this station measures a south-southeastwind at 10 m above ground (the slope of the terrain faces the south-southeast).Output for the airport is also reasonable. Environment Canada records wind direction in 10degree increments, which explains why the observed values look “chunky” (see around the 360°line on August 16th) when plotted. The model does poorest in in Terrace though this shouldnot come as a surprise. Terrace, located at the confluence of the Terrace - Kitimat valley aswell as the Skeena and Kalum rivers, is known by locals as “the bowl” because of the tendencyfor wind to blow from a multitude of directions. Annual wind roses (not shown) show thatthe dominant wind direction in Terrace is south yet in the spring and summer this is not thecase, as shown by observations in figure C.4. Statistical metrics confirm these observation, aslisted in table C.4.Table C.4: Modelled - observed wind direction statistics from the WRF evaluation.Spring SummerLocation MB MAE MB MAE(deg) (deg) (deg) (deg)Terrace B.C. Access Centre 22.603 67.736 26.230 60.731CYXT Terrace Airport 14.232 58.595 14.477 42.945Kitimat Whitesail 23.946 52.158 14.716 50.113Wind direction is a difficult parameter to evaluate statistically as the analysis requires carefulattention to detail. Given the equation for NMB it does not make sense to use this metric,as the sum of observed wind directions in the denominator does not have informational value.RMSE is also not presented, as the square of the error leads to very high values which are outof place given the mostly positive results of figure C.4. Error values (model - observed) greaterthan 180° or less than -180° were processed to confine the errors between -180° and 180° . Thefollowing is an example of the method used: If the modelled wind direction was 10 degreesand the observed wind direction was 350 degrees, modelled - observed yields a difference of-340°. The correction changed this to be a +20 degree overestimation as opposed to an -340degree underestimation.As can be seen there is a clockwise bias in the wind direction data, of between 14 and 26140degrees. This is best illustrated by figure C.4 (a) showing Kitimat’s modelled vs observedwind direction. A 20 degree clockwise bias means that if the observed wind were blowing fromthe north (360°) then the modelled wind would be blowing from the north-northeast (20°),and if the wind was observed as blowing from the south (180°) then the model would outputwind from the south-southwest (200°). MAE is high at all stations. This may be considereda poor result, yet the overall picture from the plots of Kitimat and the CYXT airport suggestthe modelled wind direction output does not appear to be unreasonable.Wind SpeedThe last parameter in the evaluation, wind speed, is also not produced directly by WRF andis calculated using the combination of U and V-wind vectors at 10m. This is done with thefollowing formula:WSPD =￿(U2 + V 2) (C.2)Model and observed wind speed are presented below for the spring and summer periods infigures C.5 (a) and C.5 (b) respectively.WRF often over predicted wind speed for Terrace (especially for the summer). The windspeeds at the airport are under-predicted but the high observed wind speeds there are alsounusual given the speeds in Kitimat (closer to the water) for the same hour. The modelperforms best in Kitimat.Table C.5: Modelled - observed wind speed statistics from the WRF evaluation.Spring SummerLocation NMB MB RMSE MAE NMB MB RMSE MAE(%) (ms-1) (ms-1) (ms-1) (%) (ms-1) (ms-1) (ms-1)Terrace B.C. Access Centre 27.212 0.573 1.704 1.331 40.482 0.854 1.825 1.386CYXT Terrace Airport -27.242 -0.969 2.401 1.885 -20.802 -0.672 2.133 1.619Kitimat Whitesail -5.760 -0.167 1.614 1.270 -5.409 -0.175 1.639 1.285As can be seen in table C.5, from the perspective of the MB, RMSE and MAE, WRF performsbest in Kitimat for wind speed. A different story is told by the NMB where speeds were over141predicted in Terrace (both seasons) by up to 40.5%, and under predicted at the airport byup to 27.2%. As well, RMSE is largest at the airport suggesting that there were occasionswhere the error was quite high, and this is confirmed by both figures C.5 (a) and C.5 (b). Itis possible that the high error values at the airport are due to measurement error, given theunusually high wind speeds measured there.142Figure C.5: Observed and modelled wind speed time series at Terrace (top), CYXT (middle)and Kitimat (bottom) for: a) spring and b) summer model periods. Black lines are observedwhile yellow, orange and blue lines are modelled.143C.2 Qualitative EvaluationThe qualitative evaluation consists of a spatiotemporal evaluation of wind streamlines andtemperature evolution over a 24-hour period. The day of July 31st was selected becausethere is a diurnal cycle in both wind direction and temperature, and the evolution of thesetwo parameters is predictable based on a basic understanding of meteorology. For the windstreamlines an offshore wind along with katabatic flow down the valley walls and in tributaryvalleys is expected in the morning. This must change at some point to an onshore wind withanabatic flow up the valley walls and in tributary valleys. For temperature, it is expectedthat the amplitude of the diurnal cycle will be greater on land than over the Douglas Channel(with respect to both heating and cooling). Mountain peaks are expected to warm as wellhowever not to the extent of the valley bottom. A more sophisticated (and quantitative) spatialevaluation is possible, however given that the concern in the TKV is advection of pollutantsup and down the valley, WRF will be deemed fit for use if it is capable of reproducing the’sloshing’ that is known to occur in the TKV from the advection of pollutants up and downthe valley over the course of one day. Differential heating of land and sea surface temperaturesas well as onshore - offshore winds are two of the main ingredients of a sea breeze, which areknown to coincide with and indeed intensify air pollution episodes (Steyn, 2003). Note thata successful spatial evaluation does not mean that the model performs optimally at all timesfor all locations, just that is capable of doing so. This approach is justified given the intendeduse of the WRF dataset and the focus on the TKV which does have monitoring stations.The locations of the Kitimat, CYXT airport and Terrace monitoring stations are not includedin the figures below so as to not obscure any of the figures. For reference on the streamlinefigure, the main north - south valley is the TKV which runs from the head of the Douglaschannel (the tip of the coastline marked in green) to approximately 3/4 up the middle of theplot. Terrace is located at the end of the TKV, at the confluence of the Skeena and Kalumrivers. The Skeena river runs from the coast through Terrace and up towards the northeastcorner of the plot. The valley that runs north from Terrace is the Kalum river valley. Better144detail can be found in figure 2.3. Seven figures are included in each section spanning midnighton the 31st to midnight on the 1st in 4-hour intervals. All times are in PST.StreamlinesAs can be seen from the first figures (especially C.6 (a)), streamlines are highly organized.Over open water (SW corner) wind is from the northwest. There are concentrated streamlinesmoving up the Skeena river, around the corner and then down the TKV. On the west side ofthe valley wind is advecting up some tributary valleys while on the east side there is katabaticmovement down valley walls and south towards the Douglas channel where it eventually meetsthe open ocean of Hecate straight. This organized pattern begins to break down by 08:00 PSTand by noon PST the arrangement of streamlines has completely reversed. At this pointwind is blowing down all tributary valleys and up the TKV towards Terrace. The ability ofWRF to resolve fine-scale topographical features can is apparent from streamlines of narrowsecondary and tertiary valleys in the domain. It is clear that 1.333 km grid spacing is requiredto accomplish this level of resolution.By the latter hours of the day streamlines have again begun to organize themselves though itcomes together in the form of strong output down the Skeena river valley (see fig C.7 (b)). Inthe TKV wind is still onshore and chaotic, especially up the valley where there appears to bea large eddy in the vicinity of the airport. Flow does not completely settle down by midnight(i.e.: became as organized as midnight the day before) but has once again turned into flowdown the TKV and Douglas Channel towards Hecate Straight.Spatiotemporal Temperature EvolutionFrom a spatial perspective the temperature resolves nicely. At midnight, as expected, thewarmest parts of the plot are all over water, while the coolest areas are over the mountainranges. Cooling continues through 04:00 PST but somewhere between C.8 (b) and C.8 (c)this pattern changes due to the early sunrise in the summer. Low-lying areas of near Kitimat145Figure C.6: Wind streamlines in the Terrace - Kitimat valley on July 31st, 2010 at: a) 00:00,b) 04:00, c) 08:00 and d) 12:00 PST. Red lines are streamlines.146Figure C.7: Wind streamlines in the Terrace - Kitimat valley on July 31st, 2010 at: a) 16:00,b) 20:00 and c) 00:00 (August 1st) PST. Red lines are streamlines.147and Terrace (Terrace is only 68 m asl) begin to heat up first as reflected in C.8 (c). Notethat the entire Skeena river valley does not warm. This is reasonable given that the valleybetween Terrace and the ocean is almost entirely water. There is very little land that is notmountain jutting out of the river and so it seems reasonable that the land around Terracewould warm faster than the Skeena river valley. By noon (fig C.8 (d)), tributary valleys havebegun to warm and the land and surrounding water are in the same temperature range, whilemountainous areas are warming but are still cool.The warmest part of the day occurs around 16:00 PST, as seen in figure C.9 (a). At this timetemperatures have reached the mid twenties in Kitimat and high twenties in Terrace. Theentire valley is warm and the land is warmer than the surrounding water. Over Hecate Straighttemperatures are still cooler which is reasonable and mountaintop temperatures remain thecoolest in the domain. The spatial graduation of temperatures from mountain top to valleybottom to ocean is reasonable. Temperatures begin to cool by 20:00 PST and by midnighttemperatures have reverted to the way they were at the start of the day with the warmesttemperatures over the Douglas channel.148Figure C.8: Spatial temperature fields for the Terrace - Kitimat valley on July 31st, 2010 at:a) 00:00, b) 04:00, c) 08:00 and d) 12:00 PST.149Figure C.9: Spatial temperature fields for the Terrace - Kitimat valley on July 31st, 2010 at:a) 16:00, b) 20:00 and c) 00:00 (August 1st) PST.150C.3 WRF Evaluation Summary and ConclusionThis section is concluded by stating that the WRF output prepared by David Suita is accept-able for use in CAMx for the three research questions. This statement is made based on thefollowing summary:• The quantitative evaluation of WRF output showed that it performed optimally withpressure and temperature, moderately optimal for wind direction and wind speed, andworst for RH (depending on the metrics used one could also argue for a poor performancefor wind speed output). For temperature the model was unable to recreate the peaks andlows of the diurnal cycles, and for wind direction there was a constant clockwise bias tothe modelled output. For wind speed the model was biased very high in Terrace (by asmuch as 40%) though in Kitimat the bias was low (by 5% which works out to be -0.167ms-1). Humidity output was greatly under predicted in Terrace by 19% in the summer(NMB). Still, even for its weakest overall performers, there were times when the outputwas reasonable and, except for RH in Terrace, the model was never truly unreasonableat any location for any parameter. The errors do create a level of uncertainty, yet thisis well within the realm of acceptable uncertainty.• The qualitative spatial evaluation showed that both wind streamlines and temperatureevolve as expected over the course of space and time for the chosen day in the evaluation.Streamlines were concentrated in overnight hours and there was evidence of katabaticflow on one side of the Terrace - Kitimat valley. The offshore breeze shifted to becomean onshore breeze about the time one would expect for a day in July, and the model wasable to resolve flow around very fine-scale topographical features like tertiary valleys.The spatial and temporal evolution of temperature was as expected for a summer dayin complex coastal terrain.151Appendix DControl Case Model Evaluation for NOand NO2D.1 Nitric OxideFrom figures D.1 (a) and D.1 (b) it is clear for both periods that modelled Control Case NOfailed to recreate observed NO for either period. This lack of characterization likely explainswhy O3 mixing ratios remain elevated overnight as there was insufficient NO to titrate O3 asper reaction 1.2. Mixing ratios are so low they are barely visible in figure D.1 and insets arerequired to magnify the modelled NO time series.The poor performance of the model is mirrored in the statistics presented in table D.1. Despitea MB of only -2.2 and -1.8 ppb for spring and summer respectively, modelled concentrationswere biased low by nearly 100% (NMB, both seasons) and there was little linear relationshipbetween modelled and observed values. These results point strongly to an improper charac-terization of NOX emissions in Smithers.152Figure D.1: Observed and modelled NO time series at Smithers for: a) spring and b) summermodel periods. Black lines are observed and red lines are modelled Control Case. Insets ineach figure allow the modelled NO to be seen.Table D.1: NO model evaluation statistics for the Control case - (all hours)Spring SummerDescription n NMB MB RMSE MAE r n NMB MB RMSE MAE r(%) (ppb) (ppb) (ppb) (%) (ppb) (ppb) (ppb)Paired mixing ratios 614 -98.386 -2.183 4.614 2.183 0.020 367 -97.738 -1.763 3.633 1.766 0.132Max model / max obs 27 -99.084 -10.756 13.874 10.756 0.032 16 -98.906 -9.390 11.433 9.390 0.531Max model w/ paired obs 27 -95.690 -2.208 2.772 2.219 0.269 16 -96.596 -2.946 4.028 2.954 0.761153D.2 Nitrogen DioxideGiven the low NO mixing ratios is is not surprising that NO2 mixing ratios also did notagree with observed values. At nighttime when there is no photolysis of NO2 reaction 1.2leads to the removal of O3 (Sillman, 2012), however without NO this is not possible. (Inunpolluted environments such as Smithers it is common for NO2 to peak overnight or in theearly morning hours, as shown in figure D.2.) During the day, reactions involving HO2 andother peroxy radicals also contribute to NO2 formaiton (Crutzen, 1979), so the lack of NO2during afternoon hours cannot be wholly explained by the absence of NO.The performance of NO2 was not as poor as NO yet the statistical metrics are not encouraging.While both spring and summer RMSE were less than 5.0 ppb, modelled NO2 was biased low by91.3 and 87.1% respectively. The maximum modelled mixing ratio was negatively correlatedto the maximum observed mixing ratio, reinforcing that, since the maxima occurred duringovernight hours, emissions were improperly characterized.Table D.2: NO2 model evaluation statistics for the Control case - (all hours)Spring SummerDescription n NMB MB RMSE MAE r n NMB MB RMSE MAE r(%) (ppb) (ppb) (ppb) (%) (ppb) (ppb) (ppb)Paired mixing ratios 615 -91.330 -3.530 4.940 3.532 0.228 366 -87.130 -2.230 3.102 2.235 0.414Max model / max obs 27 -93.129 -10.113 11.095 10.113 -0.263 16 -90.599 -6.263 6.728 6.263 -0.120Max model w/ paired obs 27 -86.221 -4.669 5.685 4.669 -0.057 16 -83.229 -3.225 3.777 3.225 0.332154Figure D.2: Observed and modelled NO2 time series at Smithers for: a) spring and b) summermodel periods. Black lines are observed and red lines are modelled Control Case.D.3 Vertical layersInvestigating model output through vertical layers in CAMx was conducted to examine if andwhen pollutants leave the surface layer and enter other layers. Figure D.3 presents plots ofPBLH, wind speed as well as contours of NO2 and NO mixing ratios for the lowest eightlayers of the model from May 6th to 12th, 2010. Note that the scales for both pollutants(located on the right hand side) have similar colours but are different in magnitude. It is clearfrom this image that all NO emitted during the day is mixed into upper layers in the model155corresponding with an increased PBLH. Overnight mixing ratios of NO and NO2 are very low,and what is present does not appear to diffuse out of the model’s surface layer. This stronglyreinforces the idea that lack of overnight O3 titration is occurring due to mischaracterizationof NOX emissions as opposed to diffusion into non-surface layers.156Figure D.3: Modelled planetary boundary layer height (top) and wind speed (second from top) time series as well as vertical NO2(second from bottom) and NO (bottom) profiles at Smithers from May 6th - 12th, 2010. Only the lowest eight layers are includedin the vertical profiles. Note that while the colour profiles of the lower plots are the same the mixing ratio scale is not.157Appendix ESensitivity Test Results for NO and NO2E.1 Nitric OxideAs seen in figures E.1 (a) and E.1 (b), the scenario which seemed to replicate the generalrange of observed NO mixing ratios (especially in the spring) was the 100xm20xp case. Whilethe timing was somewhat misplaced, this high-emitting scenario was able to produce ambientmixing ratios in the 20 - 30 ppb range, which matched what was observed on some evenings(though in the summer there were two days with erroneously high mixing ratios). The 50xmscenario was also able to produce elevated NO in the spring, while the 10xm20xp and 20xpscenarios were able to produce elevated NO, particularly in the summer (note the 10xm20xpline is directly behind the 20xp line for much of both figures). The 10xm scenario, becauseovernight emissions in this case were low, was unable to produce elevated NO.No case was able to replicate the exact timing of overnight emissions, and this affected thestatistical performance of all cases (all metrics). None of the scenarios were able to producemeaningful correlation with the observed values. All biases were low except for the 100xm20xpcase which was biased very high, with a 179% high bias in the summer period. The MAE wasonly 2 ppb or less for all cases except the 100xm20xp case, indicating that on average model158output was close to observed values (which makes sense as NO mixing ratios are near zero forlarge portions of the day) though the RMSE of 3.1 - 4.5 ppb (including spring and summerand ignoring 100xm20xp) indicates that there were some misses with peak values. Despite thelarge negative biases however, all scenarios produced improved statistics compared to the Basecase, and the NO magnitudes of the 100xm20xp case match the magnitudes of the observations(with two exceptions in the summer period on April 4th and 5th). For the two modest caseswith increased point source emissions (20xp and 10xm20xp), the biases improved from -98.4and -97.7% (spring and summer respectively) to -72.1 and -56.0%. The MAE of these twocases also improved compared to the base case, decreasing from 2.1 to 1.8 ppb in the springand 1.7 to 1.5 ppb in the summer. The complete results from this case are listed in table E.1.Table E.1: NO model evaluation statistics for the sensitivity tests - (all hours)Spring SummerDescription n NMB MB RMSE MAE r n NMB MB RMSE MAE r(%) (ppb) (ppb) (ppb) (%) (ppb) (ppb) (ppb)Base 614 -98.386 -2.183 4.614 2.183 0.020 367 -97.738 -1.763 3.633 1.766 0.13210x Mobile base Point 616 -92.901 -2.081 4.576 2.084 0.065 367 -89.158 -1.608 3.538 1.634 0.189Base Mobile 20x Point 616 -72.133 -1.616 4.419 1.801 0.013 367 -51.976 -0.938 3.265 1.505 0.19610x Mobile 20x Point 616 -71.630 -1.604 4.412 1.803 0.028 367 -55.090 -0.994 3.265 1.492 0.21350x Mobile base Point 616 -63.200 -1.416 4.353 1.791 0.069 367 -36.484 -0.658 3.195 1.607 0.235100x Mobile 20x Point 616 47.734 1.069 5.602 3.116 0.093 367 179.710 3.242 6.627 4.399 0.231159Figure E.1: Observed and modelled NO time series at Smithers for: a) spring and b) summermodel periods. Black lines are observed and red lines are modelled Control Case. Green,blue, purple, orange and light blue lines are the modelled sensitivity test cases: 10xm, 20xp,10xm20xp, 50xm and 100xm20xp respectively, as defined in table 3.5.160E.2 Nitrogen DioxideNO2 is an important gas when considering the performance of these sensitivity case studiesas NO2 is the product of the titration reaction between O3 and NO as per reaction 1.2. Itis possible for models to properly model O3 and NO, but produce elevated NO2; this can bethought of getting O3 and NO right for the wrong reasons. An example of great O3 and NOresults but poor NO2 results is the 100xm20xp scenario, as shown in E.2 (a) and E.2 (b).These plots indicate that NO2 mixing ratios are many times higher than those observed forthis extreme case. The 50xm case also produced NO2 mixing ratios well above those observed,illustrating how emissions of such magnitudes are unrealistic. A visual inspection of the timeseries show that the performance of the remaining scenarios was not uniform across bothseasons, with the 10xm20xp and 20xp scenarios matching the range of observed values in thespring and the 10xm scenario matching the range of observed values in the summer.Based on the statical metrics listed in table E.2, it can be seen that the 20xp scenario has thesmallest bias (only 8.0%) in the spring and the 10xm has the smallest bias in the summer aswell as the smallest error for both seasons. These are a considerable improvement over thebase case and a marked indicator of improved fitness. Excluding the r, these scearios were thebest performers with respect to NO2. The 100xm20xp and 50xm scenarios produced extremebias and errors, while the 10xm20xp case did not perform as well as the other modest cases.Table E.2: NO2 model evaluation statistics for the sensitivity tests - (all hours)Spring SummerDescription n NMB MB RMSE MAE r n NMB MB RMSE MAE r(%) (ppb) (ppb) (ppb) (%) (ppb) (ppb) (ppb)Base 615 -91.330 -3.530 4.940 3.532 0.228 366 -87.130 -2.230 3.102 2.235 0.41410x Mobile base Point 615 -57.243 -2.213 4.056 2.550 0.259 366 -36.669 -0.939 2.222 1.492 0.455Base Mobile 20x Point 615 7.986 0.309 4.592 3.011 0.159 366 66.170 1.694 4.567 3.001 0.35610x Mobile 20x Point 615 17.822 0.689 4.637 3.041 0.189 366 80.488 2.060 4.563 3.036 0.38650x Mobile base Point 615 101.928 3.940 7.096 5.150 0.260 366 202.092 5.173 7.256 5.325 0.467100x Mobile 20x Point 615 348.699 13.478 17.468 13.755 0.270 366 562.388 14.396 17.861 14.396 0.479161Figure E.2: Observed and modelled NO2 time series at Smithers for: a) spring and b) summermodel periods. Black lines are observed and red lines are modelled Control Case. Green,blue, purple, orange and light blue lines are the modelled sensitivity test cases: 10xm, 20xp,10xm20xp, 50xm and 100xm20xp respectively, as defined in table 3.5.162E.3 ScatterplotsThe scatterplots shown in figures E.3 and E.4 depict hourly paired mixing ratios of observedvs. modelled O3, NO and NO2 for both the spring and summer cases. The 50xm scenariowas omitted to allow the remaining 5 sets of model output to fit on one line. A linear fit isalso shown. For the O3 scenarios, daytime paired values are included and are indicated by ablack dot in the middle of the coloured circle. The strongest model performer, as indicatedby the results of the previous sections and section 3.3.2, is outlined with a blue square. Asis apparent, the 100xm20xp scenario performs well for O3 and NO but poorly for NO2. The10xm, 20xp and 10xm20xp scenarios appear to be an improvement over the base case withrespect to paired values for all gaseous pollutants (both seasons). For O3, it is clear, especiallyin the spring, that increasing emissions improves the model output, placing more and morepaired values within the 1:2 ratio lines. Daytime mixing ratios are mostly within these limits.The NO scatterplots indicate that no model scenario correctly captured emissions in Smithers,though the 100xm20xp scenario did result in elevated NO mixing ratios. The NO2 plotsindicate that the best performers were the 20xp scenario for the spring and 10xm scenario forthe summer. In order to avoid over-titration, it is proposed that model fitness is sufficientlyimproved with the inclusion of additional low-magnitude, area-based NOX increases, such asthose from the 20xp case.163Figure E.3: Observed vs. modelled O3, NO and NO2 scatterplots at Smithers for the spring model period. a), b), c), d) and e) showobserved vs. modelled O3 (Control and sensitivity cases: 10xm, 20xp, 10xm20xp, and 100xm20xp respectively). f), g), h), i) and j)show observed vs. modelled NO (Control and sensitivity cases: 10xm, 20xp, 10xm20xp, and 100xm20xp respectively). k), l), m),n) and o) show observed vs. modelled NO2 (Control and sensitivity cases: 10xm, 20xp, 10xm20xp, and 100xm20xp respectively).Hollow dots (e.g.: ) indicate all paired values while filled dots (e.g.: ) indicate daytime hours (only for O3 figures). A linear bestfit is shown as a solid black line while a dashed line represents the daytime best fit (O3 row only). 2:1, 1:1 and 1:2 ratios are shownas thin dotted lines. The best performing sensitivity case (based on a combination of all statistical metrics) is outlined in blue.164Figure E.4: Observed vs. modelled O3, NO and NO2 scatterplots at Smithers for the summer model period. a), b), c), d) and e)show observed vs. modelled O3 (Control and sensitivity cases: 10xm, 20xp, 10xm20xp, and 100xm20xp respectively). f), g), h), i)and j) show observed vs. modelled NO (Control and sensitivity cases: 10xm, 20xp, 10xm20xp, and 100xm20xp respectively). k), l),m), n) and o) show observed vs. modelled NO2 (Control and sensitivity cases: 10xm, 20xp, 10xm20xp, and 100xm20xp respectively).Hollow dots (e.g.: ) indicate all paired values while filled dots (e.g.: ) indicate daytime hours (only for O3 figures). A linear bestfit is shown as a solid black line while a dashed line represents the daytime best fit (O3 row only). 2:1, 1:1 and 1:2 ratios are shownas thin dotted lines. The best performing sensitivity case (based on a combination of all statistical metrics) is outlined in blue.165