UBC ATSC 595 ADM Homework

Homework - HYSPLIT

[ Home | Homework | ]

Homework Assignments

HW	Due by	Details
	19 Oct	Bring your laptops to class. Tim Chui will lead us in the installation of the hysplit model. If he finishes quick enough, then Reagan McKinney will start leading us thru the hysplit tutorials.
	21 Oct	Bring your laptops to class. Reagan McKinney leads us in the hysplit tutorials.
.	26 Oct	Readings: Start by viewing the first item linked in the Resources section of the Hysplit web page for this course. Next, Read the two items listed in the Theory section of the Hysplit web page for this course. Then read the BAMS journal paper about Hysplit, which is the second link in the Resources section of the Hysplit web page for this course. Also, read the Supplement to the BAMS paper. Assignment: Present to the class TODAY a section of the BAMS Supplement, as was selected during a previous class meeting. Possible sections include: Mixing Options, first paragraph (including eqs. ES7-ES8) Mixing Options, second paragraph (including eqs. ES9 - ES17) Mixing Options, third paragraph (including eqs. ES18-ES25) Stability & Mixing Depth Calculation Options Dynamic Management...for Particles and Puffs. Time-varying Emissions (first half) Time-varying Emissions (last portion, on Plume Rise)	Be prepared to discuss how mean wind and turbulence are included in Lagrangian models for atmospheric dispersion. Be prepared to discuss how Lagrangian models for dispersion in turbulent atmospheres differ from the Langevin model for molecular diffusion. Be prepared to discuss the boundary layer and the dispersion parameterizations in the BAMS supplement.
.	28 Oct	# Similarity PBL parameterizations for sigma used in hysplit # R. Stull, 23 Feb 2016 # ES = equation number in BAMS supplementary material by Stein et al (Dec 2015) on Hysplit model # http://journals.ametsoc.org/doi/suppl/10.1175/BAMS-D-14-00110.1/suppl_file/10.1175_bams-d-14-00110.2.pdf # Givens: # general tlhoriz = 10800. # lagrangian timescale for horizontal turbulence (s) g_by_tv = 0.0333 # g/Tv approx constant (m s^-2 K^-1) rho_cp = 1231. # air density times specific heat at const press, approx as constant. (W/m2)/(K m/s) # stable and neutral pbl ustar_n = 0.2 # friction velocity (m/s) for neutral PBL with moderate mean wind zi_n = 500. # pbl depth in neutral and stable tlvert_s = 5. # lagrangian time scale for vertical motions for stable pbl (s) tlvert_n = 200. # lagrangian time scale for vertical motions for unstable pbl (s) # unstable pbl Hsfc = 280. # surface heat flux (W/m2) FHsfc = Hsfc / rho_cp # kinematic heat flux (K m/s) at the surface ustar_u = 0.05 # friction velocity (m/s) for convective PBL with light winds zi_u = 1000. # mixed-layer depth for unstable pbl (m) tlvert_u = 200. # lagrangian time scale for vertical motions for unstable pbl (s) # Hint 1: for unstable pbl, use w* = Deardorff velocity from Stull Practical Meteor, eq (18.19a) # Hint 2: kinematic heat flux FHsfc = (traditional heat flux) / rho_cp . # Hint 3: for the Obukhov length (L), use z/L = -0.4 * (z/zi) * (w/u)^3 # Hint 4: the "surface layer" is generally the bottom 10% of the whole boundary layer # Exercise: # Calculate and Plot all 3 velocity variances (sigma_u2, sigma_v2, sigma_w2) in (m2/s2) # vs. height (z) for z ≤ zi, for # (1) stable/neutral pbl on one graph, using eqs (ES9 - ES11) # (2) unstable pbl on another graph, using eqs. (ES15 - ES16) # (3) unstable surface layer in the previous graph, using eq. (ES17). # Check your answers and discuss the significance. # Turn in your code and the contour plots.	This is a chance to learn about some of the planetary boundary layer (pbl) parameterizations used in Hysplit. You can write your code in any language; e.g., R, MatLab, python, fortran, excel, etc.
.	28 Oct - continued	Presentation by students on outcome of the AERMOD homework model run. Stull and/or Tim covers Panoply installation and use. Topics for in-class exercises: Get NetCDF files and view them using NASA's Panoply program. Convert WRF meteorology NetCDF files into ARL format that Hysplit can read. Then use those ARL files to run your own simulated emissions from UBC. In class, download to your laptop the zipped WRF output files , and unzip it to give a folder of NetCDF files. Add the suffix .nc to each of those filenames. (This identifies them as NetCDF files.) From the NASA Panoply website ( https://www.giss.nasa.gov/tools/panoply/ ) Mac users must use the Safari browser (not Chrome) to: first download to your laptop the updated version of Java runtime environment, follow the instructions to install Java, and to use Safari/preferences/websites to activate (using the check box) the Java plug-in, and to allow the giss.nasa.gov website to use this plug-in. Finally download the panoply application to your applications folder. Use Panoply to view the contents of any one of the NetCDF files and to view maps of different fields at different heights. This is to get you familiar with the type of info that is in NetCDF files. (Hint: to zoom any map, on the Mac use command drag to create an outline of the region you want to zoom.) Convert the WRF meteorological files from NetCDF to ARL format. If you are unable to do that conversion, then use an already converted ARL file (for a different date) provided in the next item below. If you are unable to convert the NetCDF files to ARL files yourself, then instead you can download the large (1.5 GB) ARL meteorological file we provide for Canada (at 12 km grid spacing) to your laptop, to calculate trajectories and concentrations for emissions from a hypothetical source as specified in class. (Note: it might take a few minutes to download the large data file.) Use the Hysplit gui tools to open and view the weather patterns in the ARL file. Create emission from UBC at different heights to track using the trajectory option. Also try emissions at one height to show concentrations. Also try to view a cloud of particles emitted from one source. If time, Start discussing how Calpuff works. Stull: Intro to Calpuff "slugs". Review of sections 2.1 and 2.2 of Calpuff Users Guide v5.	Here we use real weather data and learn how to set up our own emissions.
.	2 Nov	Tim Chui will show us how to install CAPUFF today in class. Meanwhile, before today, finish the following homework to turn in today. # Homework, Write your own 1-D Lagrangian particle dispersion model # R. Stull, 22 Feb 2016 # ES = equation number in BAMS supplementary material by Stein et al (Dec 2015) on Hysplit model # http://journals.ametsoc.org/doi/suppl/10.1175/BAMS-D-14-00110.1/suppl_file/10.1175_bams-d-14-00110.2.pdf # Givens: # general g_by_tv = 0.0333 # g/Tv approx constant (m s^-2 K^-1) rho_cp = 1231. # air density times specific heat at const press, approx as constant. (W/m2)/(K m/s) # unstable pbl Hsfc = 280. # surface heat flux (W/m2) FHsfc = Hsfc / rho_cp # kinematic heat flux (K m/s) at the surface zi = 1000. # mixed-layer depth for unstable pbl (m) tlvert = 200. # lagrangian time scale for vertical motions for unstable pbl (s) # plume info nparticles = 20 # number of particles (emitted as a single burst) zs = 100. # source emission height (m) for all particles # initial conditions t = 0. # initial time it = 0. # initial time index zp = zs # initial particle height for any one particle (same for all particles) wp = 0. # initial W' vertical velocity value for that one particle (same for all particles) # model integration delt = 5. # timestep increment (s) for each time step itmax = 200 # max number of timesteps # Hint 1: for unstable pbl, use w* = Deardorff velocity from Stull Practical Meteor, eq (18.19a) # Hint 2: kinematic heat flux FHsfc = (traditional heat flux) / rho_cp . # Hint 3: include particle reflection at the top and bottom boundaries. # Homework, create a Lagrangian model to track the vertical movement of a bunch particles released # simultaneously from the same source height into a convective boundary layer. Consider only vertical # particle movement, and plot the vertical position of all the particles as a function of timestep. # Use a vertical domain of 0 ≤ z ≤ zi, where zi is the mixed layer depth. # No mean wind, thus no horizontal domain. Instead, the horizontal axis in your plot is timestep index. # # Use the sigma_w expression given by Stein el al (2015) eq. ES15 in their supplement, to find the # the value of sigma_w at the particle location. # Use the position forecast eq (simplified Draxler eq 53). Z(t+∆t) = Z(t) + W' * ∆t # Make two different experiments, with slightly different eqs for forecasting W' : # # 1) Here is a simplified version of Stein et al Supplement eq. (ES4), but without the # sigma-gradient term (with notation the same as in that journal paper): # wp = (Rwp) + (sigma_w lambda * sqrt(1 - R^2) ) # # 2) Same as (1), but with the additional sigma_w gradient term dsigma_w/dz included: # wp = (Rwp) + (sigma_w lambda * sqrt(1 - R^2) ) + (sigma_w * tlvert * (1-R) * dsigma_w/dz ) # # 3) Discuss the difference between your plots from (1) and (2) above. Why is it important? # # where # tlvert = Lagrangian time scale for vertical motion, # R = autocorrelation coefficient, # lambda is a random number from a Gaussian distribution with mean of zero and st.dev of 1 # # Hint 4: Caution calculating the vertical sigma_w gradient at the top & bottom domain boundaries. # Check your answers and discuss the significance. # Turn in your code and the contour plots.	This is a chance to learn about Lagrangian particle modeling similar to what is used in Hysplit, except in only one dimension (vertical). You can write your code in any language; e.g., R, MatLab, python, fortran, excel, etc.
	4 Nov	If needed, Tim Chui will finish showing us how to install Calpuff, Calmet, Calpost, and other components. Before class, please skim the Calpuff User Instructions. Meanwhile, start working on the following homework, due 9 Nov : Homework: 2-D (x-z) Lagrangian particle transport. Use any language: R, Matlab, python, etc. Hint: You might be able to re-use parts of the code from your 1-D Lagrangian homework. ========= A) Implement the Hysplit mean-wind advection scheme in your own code, for just the U wind. 1) First: Implement the following artificial mean wind environment where the mean U ind varies in x and z, within which you may test your 1-D (x-direction) mean advection by the U component of wind. Assume V = 0 , W = 0. And neglect turbulence in the x-direction (i.e., u' = 0). I do this by assuming U obeys the log wind profile with z, and that u* varies with x. # Given: ustar0 = 0.05 # first u* constant (m/s) ustar1 = 1.0 # second u* constant (m/s) zo = 0.02 # roughness length (m) L = 50000. # horizontal domain (m). 0 ≤ L ≤ 50 km zi = 1000. # your veritical domain is 0 ≤ z ≤ zi Hint: create functions for the following equations: Use u* as a function of x: ustar = ustar0 + ustar1( (x/L) (1. - (x/L)) ) Use U as a function of z and u: U = (ustar/0.4) log( (z+zo)/zo ) 2) To see what you produced, calculate U(x,z) for fairly high resolution in x, and z, and then plot it as a contour plot. 3) Implement the Hypslit mean advection scheme, and test it for a particle released at (x, z) = (0, 100. m). Namely the source height is zs = 100. This scheme is given by eq. (1) in the Stein et al paper (2015) in BAMS. (See our course web page.) Namely, neglecting vertical turbulence for now, calculate the postion of the one particle at each time step (pick reasonable time steps) as it advects to the East, and plot the positions of the particle within the (x, z) domain. =========== B) Add the vertical turbulent dispersion code from your previous 1-D Lagrangian homework. Use your code that included reflection from z=0 and z=zi, and which prevented accumulation of particles near the top and bottom edges. # Givens: # general g_by_tv = 0.0333 # g/Tv approx constant (m s^-2 K^-1) rho_cp = 1231. # air density x specific heat at const press, is constant. (W/m2)/(K m/s) # unstable pbl Hsfc = 50. # surface heat flux (W/m2) FHsfc = Hsfc / rho_cp # kinematic heat flux (K m/s) at the surface wstar = (g_by_tv * FHsfc * zi)^(1/3) # Deardorff convective velocity tlvert = 200. # lagrangian time scale for vertical motions for unstable pbl (s) dt = 60. # timestep increment (s) for tracking any one particle # plume info nparticles = 1 # number of particles (emitted as a single burst at t = 0) zs = 100. # source emission height (m) for all particles 1) Plot the track of this 1 particle in your (x, z) domain, until it exits the domain. Use light gray colour for all the plotted points. 2) Re-do (1), but with the following special colour scheme. First 10 positions use gray. Next 2 positions, use red. Next 10 positions gray. Next 2 positions, use orange. Next 10 gray. Next 2 gold. next 10 gray. next 2 green. next 10 gray. next 2 blue. Next 10 gray. Next 2 violet. Repeat the colour sequence (grays and rainbow colours) as needed to cover the whole track of the particle until it exits the domain. ======== C) Particles and concentrations. 1) Re-do B(2), but for 2 particles released simultaneously (at t = 0, z = zs). Plot their tracks in the (x,z) domain. Also, use the identical colour scheme from B(2) for each particle, except plot two points as BLACK when the first of these 2 points is t = 120 time steps from their release time. 2) Re-do C(1), but for 50 particles released simultaneously (at t = 0, z = zs). Plot their tracks in the (x,z) domain. 3) By hand, on your graph from C(2), pick a fixed domain (i.e., draw a rectange by hand) that is 5 km wide and 200 m tall, with the base of the rectangle at z = 400 m. Position this rectangle in x such that it encloses the max number of black points with the fixed height range . Count the number of black particles in this domain. 4) Re-run C(2), which should give slightly different results because of the random numbers used for the vertical dispersion. For the exact same rectangle location as in C(3), count the number of black points. 5) Repeat (re-run C(3)) three more times, and count the number of black points for each rerun. At this point, you should have a total of 5 separate counts of black particles in that one location. ======== D) Synthesis and Discussion 1) Can each count of black points represent instantaneous concentrations of pollutants? If so, what is needed to enable that representation? 2) For the repeated realizations from part (C), do they represent an ensemble, or is each realization similar to a subsequent time for the case where each run represents a different emission time (as from a continuous emission from a point source). 3) Based on your particle counts in the rectangle for each realization, what is needed to enable those counts to represent an average concentration. (Recall that air-quality standards are based on average concentrations.) === end ===
	9 Nov	The 2-D Lagrangian HW desribed above is due today.
		(switch to the Calpuff Homework web page to see subsequent calendar dates)




7	18 Feb	Readings: Start by viewing the first item linked in the Resources section of the Hysplit web page for this course. Next, Read the two items listed in the Theory section of the Hysplit web page for this course. Then read the BAMS journal paper about Hysplit, which is the second link in the Resources section of the Hysplit web page for this course. Also, read the Supplement to the BAMS paper. Modeling (see info on our Hysplit web page): Register with NOAA ARL to gain access to the Hysplit model. See the instructions on Hysplit web page to download the tcl software, used as a graphic user interface for Hysplit on the Mac. Follow our TA's instructions on our Hysplit web page to download a pre-compiled version of Hysplit for the Mac. If you have problems, perhaps you classmates or our TA can help. Run the Hysplit4 , Trajectory , QuickStart , Run Example, as outlined in the Tips by our TA. If you did not get the same output as in the Tips, consult our TA. Run the Hysplit4 , Concentration, QuickStart, Run Example, as outlined in the Tips by our TA. If you did not get the same output as in the Tips, consult our TA. From the online Hysplit tutorial ( https://ready.arl.noaa.gov/documents/Tutorial/html/index.html ), read sections 2.1 thru 2.4 (and view the "animation" linked from section 2.4) .	Be prepared to discuss how mean wind and turbulence are included in Lagrangian models for atmospheric dispersion. Be prepared to discuss how Lagrangian models for dispersion in turbulent atmospheres differ from the Langevin model for molecular diffusion. Be prepared to discuss the boundary layer and the dispersion parameterizations in the BAMS supplement.
8	23 Feb	# Similarity PBL parameterizations for sigma used in hysplit # R. Stull, 23 Feb 2016 # ES = equation number in BAMS supplementary material by Stein et al (Dec 2015) on Hysplit model # http://journals.ametsoc.org/doi/suppl/10.1175/BAMS-D-14-00110.1/suppl_file/10.1175_bams-d-14-00110.2.pdf # Givens: # general tlhoriz = 10800. # lagrangian timescale for horizontal turbulence (s) g_by_tv = 0.0333 # g/Tv approx constant (m s^-2 K^-1) rho_cp = 1231. # air density times specific heat at const press, approx as constant. (W/m2)/(K m/s) # stable and neutral pbl ustar_n = 0.2 # friction velocity (m/s) for neutral PBL with moderate mean wind zi_n = 500. # pbl depth in neutral and stable tlvert_s = 5. # lagrangian time scale for vertical motions for stable pbl (s) tlvert_n = 200. # lagrangian time scale for vertical motions for unstable pbl (s) # unstable pbl Hsfc = 280. # surface heat flux (W/m2) FHsfc = Hsfc / rho_cp # kinematic heat flux (K m/s) at the surface ustar_u = 0.05 # friction velocity (m/s) for convective PBL with light winds zi_u = 1000. # mixed-layer depth for unstable pbl (m) tlvert_u = 200. # lagrangian time scale for vertical motions for unstable pbl (s) # Hint 1: for unstable pbl, use w* = Deardorff velocity from Stull Practical Meteor, eq (18.19a) # Hint 2: kinematic heat flux FHsfc = (traditional heat flux) / rho_cp . # Hint 3: for the Obukhov length (L), use z/L = -0.4 * (z/zi) * (w/u)^3 # Hint 4: the "surface layer" is generally the bottom 10% of the whole boundary layer # Exercise: # Calculate and Plot all 3 velocity variances (sigma_u2, sigma_v2, sigma_w2) in (m2/s2) # vs. height (z) for z ≤ zi, for # (1) stable/neutral pbl on one graph, using eqs (ES9 - ES11) # (2) unstable pbl on another graph, using eqs. (ES15 - ES16) # (3) unstable surface layer in the previous graph, using eq. (ES17). # Check your answers and discuss the significance. # Turn in your code and the contour plots.	This is a chance to learn about some of the planetary boundary layer (pbl) parameterizations used in Hysplit. You can write your code in any language; e.g., R, MatLab, python, fortran, excel, etc.
9	25 Feb	Hysplit Modeling: On your own computer, navigate to your hysplit4 directory, and create a new subdirectory called "zcaptex". Into zcaptex, copy the following meteorology binary files from the Hysplit tutorial lesson 3.1 . North American Regional Reanalysis captex2_narr.bin Weather Research Forecast Model (27 km) captex2_wrf27.bin Weather Research Forecast Model (9 km) captex2_wrf09.bin Weather Research Forecast Model (3 km) captex2_wrf03.bin NCAR Global Reanalysis Extract-1 captex2_gblr.bin NCAR Global Reanalysis Extract-2 RP198309.gbl ECMWF ERA-40 Reanalysis captex2_era40.bin In Lesson 3.1, item 1, follow the tutorial instructions click on the Set File Name button to the first file above, then click on the Create Map button to view the domain. Repeat for all the files above to view their domains. In Lesson 3.1, item 2, follow the tutorial instructions to see details about each of those files. If the resulting Simulation Log window is blank, just click on it or scroll in it to make the info appear. Confirm that the first dozen or so lines in this file give info that agrees with the summary in the tutorial. In Hysplit4, create a new subdirectory "zreanalysis". Do tutorial Lesson 3.2, item 1, by accessing Meteorology / ARL Data ftp / Reanalysis, to download that file. It will probably end up in your Hysplit4/working directory, after 20 minutes or so for the download. (The Download progress bar is flakey on the Mac, but if you click on it and then click off of it, the progress temporarily appears. Repeat on and off clicks to make it re-appear.) When the download is finished, you can move it into the zreanalysis subdirectory for future use. In Hysplit4, create a new subdirectory "zwfrt". Follow our TA's Tips to Download WRF from our own ensemble forecast system. There will be 2 netCDF files to download, both for our WRF(ARW) GFS runs: 36 km (7.5 days), and 1.3 km (3.5 days). This relates to tutorial lesson 3.3 on converting data. Pick a day for which the initial weather has boundary-layer and mid-level winds predominantly from the west (anywhere between southwest to northwest), so that we can track emissions from UBC as they move through the domain. After acquiring these netCDF files, follow the directions in tutorial lesson 3.3 to convert them to the ARL format that Hysplit can read as input. We will use these later. Read only: tutorial lessons 3.4 - 3.6. We will come back to hysplit later. Now it is time to write your own code to explore some of the parameterizations used in hysplit (see next section). .	You will gradually learn more and more about Hysplit, by working your way through their online tutorial, and using their "canned" weather and emission data.
10	1 Mar	# Homework, Write your own 1-D Lagrangian particle dispersion model # R. Stull, 22 Feb 2016 # ES = equation number in BAMS supplementary material by Stein et al (Dec 2015) on Hysplit model # http://journals.ametsoc.org/doi/suppl/10.1175/BAMS-D-14-00110.1/suppl_file/10.1175_bams-d-14-00110.2.pdf # Givens: # general g_by_tv = 0.0333 # g/Tv approx constant (m s^-2 K^-1) rho_cp = 1231. # air density times specific heat at const press, approx as constant. (W/m2)/(K m/s) # unstable pbl Hsfc = 280. # surface heat flux (W/m2) FHsfc = Hsfc / rho_cp # kinematic heat flux (K m/s) at the surface zi = 1000. # mixed-layer depth for unstable pbl (m) tlvert = 200. # lagrangian time scale for vertical motions for unstable pbl (s) # plume info nparticles = 20 # number of particles (emitted as a single burst) zs = 100. # source emission height (m) for all particles # initial conditions t = 0. # initial time it = 0. # initial time index zp = zs # initial particle height for any one particle (same for all particles) wp = 0. # initial W' vertical velocity value for that one particle (same for all particles) # model integration delt = 5. # timestep increment (s) for each time step itmax = 200 # max number of timesteps # Hint 1: for unstable pbl, use w* = Deardorff velocity from Stull Practical Meteor, eq (18.19a) # Hint 2: kinematic heat flux FHsfc = (traditional heat flux) / rho_cp . # Hint 3: include particle reflection at the top and bottom boundaries. # Homework, create a Lagrangian model to track the vertical movement of a bunch particles released # simultaneously from the same source height into a convective boundary layer. Consider only vertical # particle movement, and plot the vertical position of all the particles as a function of timestep. # Use a vertical domain of 0 ≤ z ≤ zi, where zi is the mixed layer depth. # No mean wind, thus no horizontal domain. Instead, the horizontal axis in your plot is timestep index. # # Use the sigma_w expression given by Stein el al (2015) eq. ES15 in their supplement, to find the # the value of sigma_w at the particle location. # Use the position forecast eq (simplified Draxler eq 53). Z(t+∆t) = Z(t) + W' * ∆t # Make two different experiments, with slightly different eqs for forecasting W' : # # 1) Here is a simplified version of Stein et al Supplement eq. (ES4), but without the # sigma-gradient term (with notation the same as in that journal paper): # wp = (Rwp) + (sigma_w lambda * sqrt(1 - R^2) ) # # 2) Same as (1), but with the additional sigma_w gradient term dsigma_w/dz included: # wp = (Rwp) + (sigma_w lambda * sqrt(1 - R^2) ) + (sigma_w * tlvert * (1-R) * dsigma_w/dz ) # # 3) Discuss the difference between your plots from (1) and (2) above. Why is it important? # # where # tlvert = Lagrangian time scale for vertical motion, # R = autocorrelation coefficient, # lambda is a random number from a Gaussian distribution with mean of zero and st.dev of 1 # # Hint 4: Caution calculating the vertical sigma_w gradient at the top & bottom domain boundaries. # Check your answers and discuss the significance. # Turn in your code and the contour plots.	This is a chance to learn about Lagrangian particle modeling similar to what is used in Hysplit, except in only one dimension (vertical). You can write your code in any language; e.g., R, MatLab, python, fortran, excel, etc.
11	3 Mar	Hysplit Exercise 1: Trajectories a. Execute Hysplit. (Hint, see the tips provided by our TA Matt.) b. use 3-day WRF meteorological data, already in ARL format, as provided by Matt. It is 641 MB. Copy it to your local machine for using it as input to Hysplit. You can find the meteorology data at: /scratch/data/hysplit_data/ARLDATA_MULTI.BIN Canada-Wide domain, Time_start = 16-02-29-00 Time_end = 16-03-02-23 Hint: Use the Meteorology tab in Hysplit to check if the data is good, and plot maps to see where it is. c. make a Trajectory run. (Hint: see Hysplit tutorial 4.1, but using our own data.) - use an emission start time corresponding to the first time available in the meteorology - use 9 starting locations, all at Edmonton. Assume a forest fire injects smoke at all these altitudes (m) of 10, 100, 200, 500, 1000, 2000, 3000, 4000, 5000 - run for the max hours available in the meteorology (72 hours). Forward trajectory. Default Model top 10000.0 m agl. - use default values for Vertical Motion and Output. - select the correct ARLDATA_MULTI.BIN file holding the WRF meteorological data. - Save (or Save As if you want to put the set-up info in a special file location) - Run the Trajectory model. - Display the results (use default values for all fields, except pick Meters-agl instead of Pressure) and include in the HW you turn in, for both: -- the default (most of North America) display, AND another plot -- zoomed into the closest 400 km (Hint, try setting the display to show 4 rings with 100 km radial spacing.) - Discuss the significance of the results. d. - Are there any advantages to using a lower model top? (Hint: tutorial 4.2) - Are there disadvantages or physical implications to using a lower model top? e. How does the mixed layer depth vary during the forecast? (Hint: tutorial 4.3) You can use a 3-hour increment if you wish. f. - Does our meteorological data have vertical velocities included? (Hint: tutorials 4.3 and 4.8) - If so, what is their range during the period of the forecast (i.e., the max and min values)? - What are the units, and how do these units relate to vertical speeds in (m/s)? (Hint, recall p451 of Stull, Practical Meteor.) - Does your answer make sense compared to the synoptic weather situation? (Hint, use daily weather maps from http://www.wpc.ncep.noaa.gov/dailywxmap/pdffiles.html ) g. When you look at a contour map of the surface pressure field in our meteorological input file, what do you see? h. Simulate a continuous point-source emission at 200 m height at Edmonton by emitting a particle every hour, with trajectory durations of 72 hours, for that one level. (Hint: tutorial 4.11). (I had difficulties with this.) Note: No matter what we do, Henryk has always done it first. See http://weather.eos.ubc.ca/Interpretation/StreakLines/Streaks1.CParks.mov i. Display the flow field at 200 m elevation at the start of day 1 by emitting particles from a grid of emitters. (Hint: tutorial 4.12) Hints: forward trajectories of 6 hours. Pick lower left and upper right corners for the display that are within the forecast domain (else the model crashes). j. Create an ensemble of runs from near the one emission grid point at Edmonton (200 m emission height). Show the resulting plot, and discuss the meaning of the spread of trajectories. (Hint: tutorial 4.13). (more to come soon)
12	31 Mar	Homework: 2-D (x-z) Lagrangian particle transport. Use any language: R, Matlab, python, etc. Hint: You might be able to re-use parts of the code from your 1-D Lagrangian homework. ========= A) Implement the Hysplit mean-wind advection scheme in your own code, for just the U wind. 1) First: Implement the following artificial mean wind environment where the mean U ind varies in x and z, within which you may test your 1-D (x-direction) mean advection by the U component of wind. Assume V = 0 , W = 0. And neglect turbulence in the x-direction (i.e., u' = 0). I do this by assuming U obeys the log wind profile with z, and that u* varies with x. # Given: ustar0 = 0.05 # first u* constant (m/s) ustar1 = 1.0 # second u* constant (m/s) zo = 0.02 # roughness length (m) L = 50000. # horizontal domain (m). 0 ≤ L ≤ 50 km zi = 1000. # your veritical domain is 0 ≤ z ≤ zi Hint: create functions for the following equations: Use u* as a function of x: ustar = ustar0 + ustar1( (x/L) (1. - (x/L)) ) Use U as a function of z and u: U = (ustar/0.4) log( (z+zo)/zo ) 2) To see what you produced, calculate U(x,z) for fairly high resolution in x, and z, and then plot it as a contour plot. 3) Implement the Hypslit mean advection scheme, and test it for a particle released at (x, z) = (0, 100. m). Namely the source height is zs = 100. This scheme is given by eq. (1) in the Stein et al paper (2015) in BAMS. (See our course web page.) Namely, neglecting vertical turbulence for now, calculate the postion of the one particle at each time step (pick reasonable time steps) as it advects to the East, and plot the positions of the particle within the (x, z) domain. =========== B) Add the vertical turbulent dispersion code from your previous 1-D Lagrangian homework. Use your code that included reflection from z=0 and z=zi, and which prevented accumulation of particles near the top and bottom edges. # Givens: # general g_by_tv = 0.0333 # g/Tv approx constant (m s^-2 K^-1) rho_cp = 1231. # air density x specific heat at const press, is constant. (W/m2)/(K m/s) # unstable pbl Hsfc = 50. # surface heat flux (W/m2) FHsfc = Hsfc / rho_cp # kinematic heat flux (K m/s) at the surface wstar = (g_by_tv * FHsfc * zi)^(1/3) # Deardorff convective velocity tlvert = 200. # lagrangian time scale for vertical motions for unstable pbl (s) dt = 60. # timestep increment (s) for tracking any one particle # plume info nparticles = 1 # number of particles (emitted as a single burst at t = 0) zs = 100. # source emission height (m) for all particles 1) Plot the track of this 1 particle in your (x, z) domain, until it exits the domain. Use light gray colour for all the plotted points. 2) Re-do (1), but with the following special colour scheme. First 10 positions use gray. Next 2 positions, use red. Next 10 positions gray. Next 2 positions, use orange. Next 10 gray. Next 2 gold. next 10 gray. next 2 green. next 10 gray. next 2 blue. Next 10 gray. Next 2 violet. Repeat the colour sequence (grays and rainbow colours) as needed to cover the whole track of the particle until it exits the domain. ======== C) Particles and concentrations. 1) Re-do B(2), but for 2 particles released simultaneously (at t = 0, z = zs). Plot their tracks in the (x,z) domain. Also, use the identical colour scheme from B(2) for each particle, except plot two points as BLACK when the first of these 2 points is t = 120 time steps from their release time. 2) Re-do C(1), but for 50 particles released simultaneously (at t = 0, z = zs). Plot their tracks in the (x,z) domain. 3) By hand, on your graph from C(2), pick a fixed domain (i.e., draw a rectange by hand) that is 5 km wide and 200 m tall, with the base of the rectangle at z = 400 m. Position this rectangle in x such that it encloses the max number of black points with the fixed height range . Count the number of black particles in this domain. 4) Re-run C(2), which should give slightly different results because of the random numbers used for the vertical dispersion. For the exact same rectangle location as in C(3), count the number of black points. 5) Repeat (re-run C(3)) three more times, and count the number of black points for each rerun. At this point, you should have a total of 5 separate counts of black particles in that one location. ======== D) Synthesis and Discussion 1) Can each count of black points represent instantaneous concentrations of pollutants? If so, what is needed to enable that representation? 2) For the repeated realizations from part (C), do they represent an ensemble, or is each realization similar to a subsequent time for the case where each run represents a different emission time (as from a continuous emission from a point source). 3) Based on your particle counts in the rectangle for each realization, what is needed to enable those counts to represent an average concentration. (Recall that air-quality standards are based on average concentrations.) === end ===
13	15 Apr	Hysplit Exercise 2: Concentrations ===== a. Execute Hysplit on my iMac: - run iTerm -navigate to rstull/Hysplit4/working/ -type: ./hysplit4.tcl ===== b. use 3-day WRF meteorological data, already in ARL format, as provided by Matt. It is 641 MB. If you have already copied this file to your iMac for your trajectory runs, then no need to do it again for concentration runs. Otherwise ... Copy it to your local machine for using it as input to Hysplit. You can find the meteorology data at: /scratch/data/hysplit_data/ARLDATA_MULTI.BIN Canada-Wide domain, Time_start = 16-02-29-00 Time_end = 16-03-02-23 ===== c. Do tutorial exercise 6.1, but using our own ARLDATA_MULTI.BIN data instead. (1)Menu / Reset Concentration / Setup Set start time to 16 02 29 00 1 starting location at Edmonton Run duration = 12 hours Clear old meteorology files, then add our ARLDATA_MULTI.BIN Leave all the other default values (eg. emission height is 10 m) Save (2)Concentration/Run Model Exit the running pop-up menu. (3) Concentration / Display / Concentration / Contours Execute Display Show a plot of the results, and Discuss what it shows you, and how it determines concentration from the individual particles. ===== d) Do this run again, but for 36 hours. Plot the concentrations. Compare with the plot from (c). to answer these questions: - what is the emission duration - what is the averaging period - why does the last plot look like particles instead of concentrations? ===== e) Do tutorial exercise 6.2 (parts 1-4), but using our own ARLDATA_MULTI.BIN data instead. (i) Use the same set-up as (c1), namely, use 12 h run and 10 m emission. Plot the result and discuss. (ii) Repeat, but with 36 h run. Compare the results. Discuss what you see. Also, based on what you now know, you might want to change your discussion to questions (c) or (d) above. (iii) Repeat part (i), but only for 12 h, and for the Grids set up for an averaging layer depth of 1000 m. Discuss why the concentrations values are different from (3(i)). ===== f) Do Tutorial 6.3, but with our own arl file. Start with the same setup as you just finished in the previous exercise. (i) Do Tutorial 6.3 parts (1) and (2) and Discuss how this relates to what you saw in your previous exercises. (ii) For Tutorial 6.3 part (3), if you can't create an animation, don't worry. You can view the pdf output on your mac in Preview. Set it to view single frame, then use the arrow buttons on your keyboard to step forward or back. (ii) Do Tutorial 6.3 part (4). Namely: Concentration/SetupRun/SaveAs and: Advacnced/ConfigurationSetup/Concentration/SaveAs ===== g) Do Tutorial 6.4, but with our ARL file. (i) Do tutorial step (1) and in Concentration / Setup change the run time to 36 hours before you run the model as step (2). Show and discuss the results. (ii) Do tutorial step (3), and show the resulting contour plots. Do tutorial step (4) , but changing the grid spacing to 0.2 degrees in latitude and longitude (not to 0.005 as stated in the tutorial). Show the results. - Discuss how averaging time and averaging volume affect the concentration that is computed, and how it is displayed. For example, compare the concentrations (max value, and shape of distribution) at 1800 on 29 Feb from the output of tutorial (4) vs. tutorial (3). - Further discuss whether the concentration SHOULD change when the time and spatial averaging is changed? Why? - In the last run that you contoured, what is the red dot in the graphics output? - Also discuss how the apparent run time relates to the number of particles released. - Further discuss what is the "best" averaging time and averaging volume to use for forecasts of real events where the actual concentration values are important. How do you determine this "best" value? Save the control and namelist files again, with the filenames suggested at the end of tutorial 6.4.
14	29 Apr	Hysplit Exercise 3: Concentrations (continued) ======== h) Do tutorial 6.5. Setup using our own ARL meteorological data, with emissions from Edmonton at the default 10 m height, and for only a 12 hour forecast. In the Advanced Config Setup menu, for menu #4, change the particles released per cycle back to 100. Run using the SETUP.CFG file. Follow the tutorial for the display, but show every 1 particle, and zoom about 70%. For that same output, also produce the normal contoured concentration display. Show both displays, (1) Compare and Discuss the results. (2) What do the different colour particles represent in the particle display? How does the vertical cross section relate to the movement of particles in the horizontal map display? What does the diagonal red dashed line represent in the particle map display? (3) Why does the normal contoured concentration display show contours that do NOT agree with the particle positions from the particle map? (4) Recall your own Lagrangian 2-D model that you wrote. Do those runs help you interpret the output from Hysplit? If so, how? ======== i) Do tutorial 6.6. Change the run duration to 36 h. And emit 10 particles each period. For the Display, try 50% zoom but if the particles leave the displayed domain then change to 30% zoom, and plot concentration contours. Note, this method is a "hybrid" method, in that the horizontal dispersion is as a puff, but the vertical is still as particles. (1) Show the results and discuss. (2) Compare the top-hat vs the Gaussian horizontal puff. In particular, why does there appear to be a different number of puffs in the Gaussian case at 36 h into the forecast compared to the top-hat case? (3) Compare the Gaussian horizontal puff output to your particle output from homework part (d) (i.e., from tutorial 6.1), particularly at the 36 hour forecast. (4) Compare and discuss the averaging times used to create your answer to part (d) to the averaging times used for the puff forecast. (5) In tutorial question 5 where we increase the particles released to 100. When you run the model, watch the simulation log output. Why does the percent complete increase rapidly at first, but more slowly as the model nears 100% completion (Hint, come back and answer this after completing tutorial 6.7)? How do the shapes and concentrations of the overall pollutant concentration pattern compare with the ones for particles that you did earlier? Which version is correct: puffs or particles? (6) Make an additional run, this one using Gaussian horizontal puffs and Top-hat vertical puffs. This is NOT a hybrid mod, but purely uses puffs in all dimensions. How does this output compare to the other puff runs you just made? Which looks more realistic? Why do the first few frames of output show NO pollutants (Hint, come back to this after you have read tutorial 6.8) ? ========= j) Read, but do not run, tutorial 6.7. Use this info to help answer the previous question i(5). ========= k) Read tutorial 6.8. (1) How do the eqs. described here compare to what you used in your Lagrangian 2-D code that you wrote? (2) Why does the eq. for concentration in exercise 5 of tutorial 6.8 have "delta_c" instead of just "c" on the left side of the eq.? (I might have asked this before.) ============ l) Do tutorials 6.9 & 6.10 but for our ARL data, and only out to 36 hour forecast. Use the start time and location as we have already been using for Edmonton. Pretend that the Sudbury smelter is located in Edmonton. The annual emission rate of SO2 was 2 million tonnes/yr decades ago when the smelter was at its peak production. Convert this to emissions in g/s, to use in Hysplit for this exercise. Assume as stack-top height of 380 m (which was the actual height of the Sudbury stack), but located in Edmonton for this exercise. Try to follow as much of this tutorial as possible, modifying as needed to work with our location, meteorology, and emissions. Show and Discuss your results. Help your classmates. ====== m) Read tutorial 6.11. Use that info and use the output from the previous model run to plot a meteogram of SO2 concentration at the surface (i.e., at z = 10 m) vs. time (0 to 36 h) for a single receptor station located 100 km East of Edmonton smelter. Show and Discuss your results, and explain why or why not it agrees with the contour plots that you had produced in the previous exercises. === You are finished with Hysplit. Note that the tutorials cover many additional features and options that we don't have time to cover. Please skim the Tutorial index (via the Home link in any of the tutorials) to see the list of features and options. But you have sufficient background to cover them on your own, if you need them for your own research. It always helps for you to create your own homework exercises for simplified scenarios, to understand the basics of any model before you apply it to real situations. ====== -end- Thanks again to Matt, our TA, for teaching us the basics, and for creating the WRF ARL file for these exercises.

Other Info

[ Home | Homework | ]