BUOYANT GAUSSIAN PLUME MODELING OFATMOSPHERIC DISPERSION OF Cs-137 EMITTED FROM A NUCLEAR REACTOR BY JOSEPH ATTEH DJANGMAH (10362178) This thesis is submitted to the University of Ghana, Legon in partial fulfillment of the requirement for the award of MPhil Computational Nuclear Sciences and Engineering degree JULY, 2013 University of Ghana http://ugspace.ug.edu.gh ii DECLARATION I hereby declare that with the exception of references to other people‘s work which have duly been acknowledged, this compilation is the result of my own research work and no part of it has been presented for another degree in this University or elsewhere. ……………………………………. Date…………………………… DJANGMAH JOSEPH ATTEH I hereby declare that the preparation of this project was supervised in accordance with the guidelines of the supervision of Thesis work laid down by the University of Ghana. ………………………………… …………...…………………………………… DR. INNOCENT J.K. ABOH NANA (Prof.) AYENSU GYEABOUR I (PRINCIPAL SUPERVISOR) (CO-SUPERVISOR) Date………………………….. Date…………………………………. University of Ghana http://ugspace.ug.edu.gh iii ABSTRACT A two dimensional Advection-Diffusion equation was developed to model atmospheric dispersion of Cs-137 in order to visualize the phenomenon of the transport ofCs-137 in the atmosphere. The study was also aimed at analyzing the activity concentration of Cs- 137 and examined the mathematical approaches for determining the mass concentration of Cs-137 released to the atmosphere. The Buoyant Gaussian model was examined and a finite difference scheme was developed to solve the partial differential equations of radionuclide pollution transport in the atmosphere. The model program code for Gaussian Plume was implemented by using the MATLAB interface contouring software. The result outputs showed that like other parameters, depth has significant influence on the dispersion of pollutants; the highest and least activity concentrations of Cs-137 dispersion were found to be 0.4686 𝑚𝑔/𝑚3and 0.0002 𝑚𝑔/𝑚3 at 2.5 𝑚/𝑠 respectively. The maximum distance that was covered by Cesium-137 from the source downwind was 1621m at 5m/s. In the event of an accidental atmospheric release of radionuclides from a nuclear power plant, accurate forecasting of the dispersion and the activity concentrations of radionuclides are required for the preparation of adequate countermeasures and evacuation. University of Ghana http://ugspace.ug.edu.gh iv DEDICATION To you, my beloved wife—Mary Boateng Djangmah and my adoring Daughter— Mehetabel Djangmakie Djangmah, I dedicate this thesis. University of Ghana http://ugspace.ug.edu.gh v ACKNOWLEDGEMENTS First and foremost, I humbly thank the Almighty God for granting me this rare opportunity to join this course of study and for his fabulous blessings all through the study. Lord, I am most grateful. I wish to extend my profound gratitude to my supervisor, Dr. Innocent Aboh for his scholarly guidance and support. I would also like to express my sincere gratitude to my co-supervisor and the course coordinator, Nana (Prof.) A. Ayensu Gyeabour I for his scholarly and systematic guidance towards the end of this work. I also appreciate the support and guidance received from my wife, daughter and all my siblings. My dear colleagues at the University played an important role; they were cooperative during the entire period of my studies. To them all, my sincere thanks. ` University of Ghana http://ugspace.ug.edu.gh vi TABLE OF CONTENTS DECLARATION ........................................................................................................... ii ABSTRACT ................................................................................................................. iii DEDICATION .............................................................................................................. iv ACKNOWLEDGEMENTS ........................................................................................... v TABLE OF CONTENTS .............................................................................................. vi LIST OF TABLES ........................................................................................................ ix LIST OF FIGURES ....................................................................................................... x LIST OF SYMBOLS AND CONSTANTS .................................................................. xi CHAPTER ONE .......................................................................................................... 1 INTRODUCTION........................................................................................................ 1 1.1. BACKGROUND OF RESEARCH ..................................................................... 1 1.2. RESEARCH PROBLEM STATEMENT ........................................................... 3 1.3 JUSTIFICATION OF RESEARCH ..................................................................... 3 1.4. OBJECTIVES OF THE RESEARCH ................................................................ 3 1.5. DESCRIPTION OF GHARR-1 .......................................................................... 4 1.5.1. Miniature Neutron Source Reactor (MNSR) ................................................ 4 1.5.2. General Site Description ............................................................................... 6 1.5.3. Weather Conditions ...................................................................................... 6 1.6. RADIOISOTOPE ATMOSPHERIC DISPERSION ACCIDENT SCENARIOS .................................................................................................................................... 8 1.6.1. Atmospheric Dispersion of Radioactive Materials ....................................... 8 1.6.2. Beyond Design Basis Accident (BDBA)...................................................... 8 1.6.3. Seismic Accident .......................................................................................... 9 1.6.4. Human Error ................................................................................................. 9 1.7. CHOICE OF Cs-137 FOR THE STUDY ........................................................... 9 1.8. ORGANIZATION OF THESIS ........................................................................ 10 CHAPTER TWO ....................................................................................................... 11 LITERATURE REVIEW ......................................................................................... 11 2.1. GENERAL DISPERSION PLUME MODELS ................................................ 11 2.1.1 Box Dispersion Model ................................................................................. 11 2.1.2. Lagrangian Dispersion Model .................................................................... 12 2.1.3. Eulerian Dispersion Model ......................................................................... 13 University of Ghana http://ugspace.ug.edu.gh vii 2.1.4. Gaussian Dispersion Model ........................................................................ 14 2.2. SYSTEM OF RADIOLOGICAL PROTECTION ............................................ 15 2.3. NUMERICAL PLUME MODELS ................................................................... 16 2.4. ADVECTION-DIFFUSION EQUATION SOLUTION METHODS ........... 18 2.4.1. Diffusion Equation Solution Methods ........................................................ 18 2.5. AIR DISPERSION MODELING OF POLLUTANTS IN THE ATMOSPHERE ....................................................................................................... 21 2.5.1. Turbulence .................................................................................................. 22 2.5.2. Mechanical Turbulence .............................................................................. 22 2.5.3. Thermal Turbulence ................................................................................... 22 2.6. EFFECTS OF WIND FLUCTUATIONS ON DISPERSION .......................... 23 2.7. FACTORS AFFECTING ATMOSPHERIC DISPERSION OF POLLUTANTS .................................................................................................................................. 23 2.7.1. Atmospheric Stability ................................................................................. 23 2.7.2. Weather Conditions .................................................................................... 24 2.7.3. Topography ................................................................................................. 26 2.8. DISPERSION PARAMETERS .................................................................... 27 2.9. FINITE DIFFERENCE METHOD FOR SOLVING THE ADVECTION- DIFFUSION EQUATION (ADE) ........................................................................... 28 CHAPTER THREE ................................................................................................... 30 RESEARCH METHODOLOGY ............................................................................. 30 3.1. PHYSICAL DESCRIPTION OF THE TRANSPORT PHENOMENON ........ 30 3.2. DERIVATION OF THE ADVECTION-DIFFUSION EQUATION ............... 32 3.3. ANALYTICAL SOLUTION OF THE ADVECTION-DIFFUSION EQUATION ............................................................................................................. 34 3.3.1. Assumptions and Boundary Conditions [38]; ............................................ 34 3.3.2. REDUCTION OF ADE TO BOUNDARY VALUE PROBLEM (BVP) TRANSFORMATION ......................................................................................... 36 3.3.3. SOLUTION OFTRANSFORMED ADVECTION-DIFFUSION EQUATION.......................................................................................................... 38 3.4. FINITE DIFFERENCE SOLUTION OF THE ADVECTION-DIFFUSION MODELFOR RADIONUCLIDE (Cs-137) ATMOSPHERIC TRANSPORT ...... 41 3.4.1. Choice of Method ....................................................................................... 42 3.4.2. Crank- Nicholson Solution of the ADE ...................................................... 42 University of Ghana http://ugspace.ug.edu.gh viii 3.4.3. NUMERICAL FLOW CHART FOR THE ADE OF THE GAUSSIAN MODEL ................................................................................................................ 45 3.4.4. STRUCTURED PROGRAM ..................................................................... 46 3.6. SOLUTION ALGORITHM BASED ON MATLAB CONTOURING ............ 47 SOFTWARE ......................................................................................................... 47 3.6.1. PROGRAM IMPLEMENTATION ............................................................ 49 3.6.2. COMPUTATIONAL ERRORS ................................................................. 49 3.6.3. CALIBRATION, VALIDATION AND VERIFICATION........................ 50 CHAPTER FOUR ...................................................................................................... 51 RESULTS AND DISCUSSIONS .............................................................................. 51 4.1. RESULTS.......................................................................................................... 51 4.1.1. SOURCE EMISSION AT GROUND–LEVEL ......................................... 52 4.1.2. CONCENTRATION OF CESIUM-137 AT ELEVATED SOURCE ........ 53 4.1.3. CESIUM-137 DISPERSION AT A SET OF RECEPTOR LOCATIONS 56 4.1.5. Analytical and Numerical Integration of Cs-137 Concentrations .............. 63 4.1.6. ERROR ANALYSIS .................................................................................. 66 4.2. DISCUSSION ................................................................................................... 67 CHAPTER FIVE ....................................................................................................... 69 CONCLUSIONS AND RECOMMENDATIONS ................................................... 69 5.1. CONCLUSIONS ............................................................................................... 69 5.2. RECOMMENDATIONS .................................................................................. 70 REFERENCES ........................................................................................................... 71 APPENDICES ............................................................................................................ 76 APPENDIX I ............................................................................................................ 76 APPENDIX III ......................................................................................................... 80 APPENDIX IV ......................................................................................................... 82 APPENDIX V .......................................................................................................... 83 University of Ghana http://ugspace.ug.edu.gh ix LIST OF TABLES 1.1. Radioactive Released in the Environment 8 2.1. Atmospheric Stability Classification Developed by Pasquill 24 2.2. Parameter Values Used in the Gaussian Plume Model 28 4.1. Mass (𝑚𝑔)of Cs-137 Deposited at each Receptor Zone:𝑢 = 2.5𝑚/𝑠 57 4.2. Mass(𝑚𝑔) of Cs-137 Deposited at each Receptor Zone: 𝑢 = 5𝑚/𝑠 59 4.3. Mass (𝑚𝑔) of Cs-137 Deposited at each Receptor Zone: 𝑢 = 7𝑚/𝑠 61 4.4. Mass (𝑚𝑔) of Cs-137 Deposited at each Receptor Zone :𝑢 = 10𝑚/𝑠 63 4.5. The relationship between the numerical and analytical solution of Cs-137 concentration level 64 4. 6. Error Analysis for ASOL & NSOL 67 University of Ghana http://ugspace.ug.edu.gh x LIST OF FIGURES 1.1. A schematic diagram of the GHARR-1 reactor vessel 5 1.2. Area Map of Kwabenya-Accra, Ghana 7 2.1. A schematic representation of the input-output of an air dispersion model 21 2.2. Wind Rose for Accra 26 3.1. Gaussian Shape of the Plume Cross-Section Centerline 31 3.2. Dispersion Control Volume of the ADE 32 3.3: Flow Chart for the ADE of the Gaussian Plume Model 45 4.1. Cs-137 Dispersion Diagram at Ground-Level 52 4.2. Cs-137 Dispersion Diagram for Elevated Height: Vertical Cross-section 54 4.3. Cs-137 Dispersion Diagram for Elevated Height: Horizontal Cross-section 55 4.4. Dispersion of Cs-137 Concentration at Multiple Receptor Zones: 56 4.5. Mass of Cs-137 (mg) Deposited at each Receptor Zone: 𝑢 = 2.5𝑚/𝑠 56 4.6. Dispersion of Cs-137 Concentration at Multiple 𝑢 = 5𝑚/𝑠 58 4.7. Mass of Cs-137 (mg) Deposited at each Receptor Zone: 𝑢 = 5𝑚/𝑠 58 4.8. Dispersion of Cs-137 Concentration at: 𝑢 = 7𝑚/𝑠 60 4.9. Mass of Cs-137 (mg) Deposited at each Receptor Zone: 𝑢 = 7𝑚/𝑠 60 4.10. Dispersion of Cs-137 Concentration at: 𝑢 = 10𝑚/𝑠 62 4.11. Mass of Cs-137 (mg) Deposited at each Receptor Zone: 𝑢 = 10𝑚/𝑠 62 University of Ghana http://ugspace.ug.edu.gh xi 4.12. A Graph of Numerical Versus Analytical Along the North-East Direction 65 LIST OF SYMBOLS AND CONSTANTS 𝜌 Particle density (mg/m3) 𝜇 Viscosity of air (mg/ms) 𝑢 Wind speed (m/s) 𝑟 Particle radius (m) 𝑔 Gravitational acceleration (m/s2) 𝑎, 𝑏 Diffusion parameters (m2) 𝑄 Emission rate (mg/s) 𝐾 Eddy diffusivity coefficient 𝐻 Stack height (m) Bq Becquerel (unit of radioactivity) 𝑉𝑠𝑒𝑡 Settling velocity (m/s) 𝑉𝑑𝑒𝑝 Deposition velocity (m/s) University of Ghana http://ugspace.ug.edu.gh xii Mathematical Symbols 𝛿 Delta 𝜕 Partial differential 𝜗 Theta variant ∞ Infinity Greek Symbols 𝜁 Zeta 𝜉 Xi 𝜂 Mu 𝜎 Sigma Variant ℒ Laplace Transform ∆ Increment University of Ghana http://ugspace.ug.edu.gh xiii ABBREVIATIONS ABS Absolute ASOL Analytical Solution BDBA Beyond Design Basis Accident DBA Design Basis Accident FDE Finite Difference Equation FDMS Finite Difference Method Solution GHARR-1 Ghana Nuclear Research Reactor 1 IAEA International Atomic Energy Agency ICRP International Commission on Radiological Protection INES International Nuclear and Radiological Event Scale MATLAB Matrix Laboratory MNSR Miniature Neutron Source Reactor NISA Nuclear and Industrial Safety Agency NNRI National Nuclear Research Institute PDE Partial Differential Equation SAR Safety Analysis Report University of Ghana http://ugspace.ug.edu.gh 1 CHAPTER ONE INTRODUCTION 1.1. BACKGROUND OF RESEARCH Nuclear reactors are built with many inherent safety features. Notwithstanding the numerous safety features promulgated by regulators in the nuclear field, accidents do occur although the occurrence is very rare. There have been three major reactor accidents in the history of civil nuclear power generation, namely the — Three Mile Island (USA, 1979), Chernobyl (Ukraine, 1986) and Fukushima (Japan, 2011) [1]. Of all the accidents and incidents, only the Chernobyl and Fukushima accidents resulted in emitting very large amount of highly fractionated radioactive debris into the atmosphere [1, 2]. The more dosimetrically significant radionuclides released to the environment were 131I, 134Cs and 137Cs [1, 2]. The determination of emission rates of various radioisotopes and prediction of pollutants concentration are necessary to assess the impacts of radionuclides on the environment. The study of the dispersion of radionuclides also determines the occurrence and frequencies of worst accident scenarios and in the end; enables regulators to minimize emissions during these adverse conditions [2]. Atmospheric dispersion modeling is one mechanism that can be used to investigate radionuclides emissions from nuclear reactors and dispersion into the environment. Atmospheric dispersion modeling is the mathematical simulation of the dispersion of pollutants in the boundary layer of the atmosphere [2, 3]. University of Ghana http://ugspace.ug.edu.gh 2 The dispersion plume models are usually based on a solution of the advection- diffusion equation by either analytical or numerical methods. Also, computer programs and algorithms can be developed to solve the mathematical equations to simulate the dispersion of pollutants [3]. A number of air pollutant dispersion models have been developed for predicting ground level concentration, in most cases for regulatory purposes. Dispersion models, however, differ in assumptions and structure as well as in the algorithms used; therefore, the predictions vary from model to model. There are several basic mathematical algorithms some of which include: Box models, Buoyant Gaussian model, Lagrangian model, Eulerian model and the Dense Gas model, but the models differ in the type of pollutant accommodated, pollutant source type and the nature of the plume been used [3].The three types of air pollution emission plumes are:  Buoyant plumes, where the plumes are lighter than air because of the higher temperature and lower density exhibited than the surrounding ambient air.  Dense gas plume, where the plumes are heavier than air, and the plumes exhibit higher density than the surrounding ambient air.  Passive or neutral plumes are plumes which are neither lighter nor heavier than the surrounding ambient air. The purpose of atmospheric transport and diffusion calculations are to provide estimates of concentration and surface deposition from routine and accidental releases of pollutants to the atmosphere. These calculations provide the link between emissions to the atmosphere and direct or indirect pathways to man for dose calculations [4]. University of Ghana http://ugspace.ug.edu.gh 3 1.2. RESEARCH PROBLEM STATEMENT The research problem addressed in this thesis was the simulation of the Buoyant Gaussian Plume to predict dispersion of radionuclides and to determine the resulting ground-level concentrations following accidental release of radionuclides into the atmosphere. 1.3 JUSTIFICATION OF RESEARCH One of the first challenges in the history of air pollution modeling was the understanding of the dispersion properties of radionuclide emitted from nuclear facilities [4], and for this purpose; the Gaussian Plume Model was developed. This model was applied for the main purpose of calculating the maximum impact from the source of the radionuclide released. The Gaussian plume model which was formulated by determining the horizontal and vertical spread of the plume, measured by the standard deviation of the plume‘s spatial concentration distribution [4].The research was therefore aimed at predicting the dispersion and ground level concentration of Cs-137 with reference to the GHARR-1 MNSR during supposed accidental situations by analytical and numerical solutions of the Advection-Diffusion Equation (ADE). Consequently, the steady-state ADE of the Gaussian plume dispersion model was solved to determine the dispersion of Cs-137 from the GHARR-1 facility in case of accidental release (and to inform Decision Makers for the preparation of adequate counter measures and possible evacuation). 1.4. OBJECTIVES OF THE RESEARCH The primary aim of the research was to adopt the Buoyant Gaussian Plume model to predict dispersion of Cs-137 and the resulting ground-level concentrations under University of Ghana http://ugspace.ug.edu.gh 4 different stability conditions for a given pollution source and the meteorological conditions that prevailed during the study. The objectives of the study were to:  develop a mathematical model of Advection-Diffusion Equation for Cs- 137 dispersion;  estimate the concentration and surface deposition from accidental release of Cs-137 to the atmosphere;  determine the mass of Cs-137 discharged at different receptor zones. 1.5. DESCRIPTION OF GHARR-1 1.5.1. Miniature Neutron Source Reactor (MNSR) GHARR-1 is a 30kW tank-in-pool type research reactor. It is cooled by natural convection and can be operated at a maximum thermal neutron flux of 1 × 1012𝑛/ 𝑐𝑚2s. The reactor complex contains five major components. These are the reactor assembly, control console, auxiliary systems, irradiation system and the pool. The reactor assembly consists of the reactor core, beryllium (Be) reflector, small fission chambers for detecting neutron fluxes, 1 central cadmium (Cd) control rod and its drive mechanism, and thermocouples for measuring inlet and outlet temperatures of the coolant. There are five inner irradiation tubes installed within the beryllium annulus. Five outer irradiation tubes are also installed outside the beryllium annulus. The reactor vessel is a cylindrical aluminium (Al) alloy container, 0.6 m in diameter and 5.6 m high. The container, which is built in 2 sections, is suspended in a stainless steel-lined water pool made of reinforced concrete [5]. The core consists of fuel elements, which form a fuel cage. The cage is inside an annular beryllium reflector and rests on a lower beryllium reflector plate. The volume of the vessel is 1.5𝑚3. The University of Ghana http://ugspace.ug.edu.gh 5 fuel elements are all enriched uranium-aluminium (U-Al) alloy extrusion clad with aluminium. The fuel elements are arranged in 10 multi-concentric circle layers at a pitch distance of 10.95 mm. The element cage consists of 2 grid plates, 4 tie rods and a guide tube for the control rod.Screws connect the 2 grid plates and 4tie rods. The total number of lattice positions is 354 and the number of fuel elements is 344. The remaining positions are occupied with 6 dummy aluminium elements as shown in Figure 1.1 [5]. Figure 1.1: A schematic diagram of the GHARR-1 reactor vessel [6] University of Ghana http://ugspace.ug.edu.gh 6 1.5.2. General Site Description The National Nuclear Research Institute (NNRI), the operating organization for GHARR-1, lies on the south-east flank of the Akwapim ridge and it is characterized by a general increase in level from south-east to north-east, though the topography as a whole is gently undulating. The broad valley of the Onyasia stream flanks the site on its eastern side; swampy conditions are in the north-east of the site, 40.0 m above sea level. Figure 1.2 shows the area map of Kwabenya where GHARR-1 is located [5]. 1.5.3. Weather Conditions The prevailing weather conditions played a major role in determining the diffusion, direction, distribution and transportation of atmospheric pollutants, hence the weather situations of an area under study must be considered. 1.5.3.1. Wind The wind speed ranged from 2.5 km/hr to 8.5 km/hr with an average wind speed of 10 km/hr, measured at 2𝑚 high [5]. 1.5.3.2. Temperature The mean minimum temperature in the area was 22℃ . The mean maximum temperature was 34℃. 1.5.3.4. Rainfall There are two rainfall seasons in Ghana. The major rainfall season is from April to July with a mean monthly rainfall total of 122 mm. The minor season is in September to October with a mean monthly rainfall total of 70 mm. The minimum annual rainfall total is 964 mm [5]. University of Ghana http://ugspace.ug.edu.gh 7 Figure 1.2: Area Map of Kwabenya, Accra-Ghana [6]. University of Ghana http://ugspace.ug.edu.gh 8 1.6. RADIOISOTOPE ATMOSPHERIC DISPERSION ACCIDENT SCENARIOS 1.6.1. Atmospheric Dispersion of Radioactive Materials Airborne radioactive released into the environment under normal operation from the GHARR-1 MNSR is mainly 42Ar from the activated air in the irradiating tubes. This is estimated to be less than 3 × 1010Bq per year. Under accident conditions, the estimated quantities of airborne radioactivity from various radionuclides released into the atmosphere are presented in Table 1.1, [5]. Table 1.1 Radioactive Released in the Environment [5] Nuclides Annual Discharge ( Bq) 90Sr 2.3 × 103 132I 6.2 × 103 137Cs 7.5 × 103 133Xe 1.1 × 105 134I 1.5 × 106 131I 1.3 × 107 131Xe 3.2 × 107 1.6.2. Beyond Design Basis Accident (BDBA) The BDBA, also referred to as the Maximum Hypothetical Accident, is not expected to occur. The BDBA is therefore not analyzed, but is described for purposes of University of Ghana http://ugspace.ug.edu.gh 9 emergency planning, and is always an accident more severe than the Design Basis Accident (DBA). The accident scenario emergency is considered with the following assumptions [5]: 1. The GHARR-1 building collapses; 2. The reactor vessel water and the pool water leak at a rate 4 m3/hr; 3. The reactor core is exposed to air after 6 hours; 4. The reactor was operating at 30 kW. 5. Aircraft accident 1.6.3. Seismic Accident The vessel is an important barrier against the release of fission products. To ensure its safety and reliability during any earthquake, the design of the vessel has taken into consideration an earthquake load. Two types of earthquakes stipulated by the United States Department of Energy (DOE) were adopted for the design. These are: Safe Shutdown Earthquake (SSE) and Operation Base Earthquake (OBE). The SSE usually takes the intensity one grade ± more than that recorded in local history while the maximum ground acceleration of OBE should not be smaller than half of the SSE [5]. 1.6.4. Human Error Human errors in critical activities such as adding beryllium reflector shim plates, spent fuel core removal and fresh fuel core installation have been considered as part of the safety analysis report [6]. 1.7. CHOICE OF Cs-137 FOR THE STUDY Although Table 1.1 refers to release fractions of some radioisotopes, for the current study, the author has decided to investigate only one radioisotope, Cs-137. The University of Ghana http://ugspace.ug.edu.gh 10 release fraction of Cs-137 is very high and volatile as compared to the rest of the radioisotopes [6]. The focus of the study is to demonstrate the dispersion of radionuclides into the atmosphere during nuclear accident, hence the reason for choosing one significant radioisotope. Therefore, Cs-137 which falls under the highly volatile group and a large contributor to the dose of an exposed person during accident has been chosen for the study. 1.8. ORGANIZATION OF THESIS Adequate background information on the underlying principle for the research is discussed in Chapter one. Also, the objectives and justification for the research are briefly discussed and a description of the GHARR-1 MNSR is highlighted in this chapter. Chapter two focuses on the review of literature pertaining to the study. The Chapter gives an overview of current literature concerning plume models, numerical solutions to the Advection-Diffusion Equation and a review of the dispersion model used in the study. In the third Chapter, all the methodologies to the governing equations and the finite difference of the Crank Nicholson approach are presented. Further, the modeling process and the research for data collection are explained. The results of the model simulation equations are presented in Chapter Four. The Chapter also provides a Discussion of the results in relation to other published work. Finally, Conclusions of the research findings are presented in Chapter Five and recommended areas to be investigated in further studies are outlined. University of Ghana http://ugspace.ug.edu.gh 11 CHAPTER TWO LITERATURE REVIEW Several studies on atmospheric pollutants dispersion using mathematical models have been carried out by various researchers all having the same goal of how atmospheric air quality can be preserved [7]. This literature review examines current research in the areas of atmospheric pollution modeling and methods of finding numerical and analytical solutions of the Advection- Diffusion Equations, (ADE). The literature reviewed in this work undertakes a broad survey of plume modeling. In this section, existing numerical plume models are discussed which also includes the research that has been done on computational modeling of one-dimensional advection-diffusion equation. 2.1. GENERAL DISPERSION PLUME MODELS The modeling of pollutants dispersion in the atmosphere is carried out by using mathematical algorithms. There are several dispersion models, some of which include: Box model, Gaussian model, Lagrangian and Eulerian model. These models differ in the type of pollutant accommodated, pollutant source type and whether they use plume or puff approach [7]. 2.1.1 Box Dispersion Model The Box model is the simplest of all modeling algorithms which is based on the conservation of mass [7]. The air shed is treated as a box into which pollutants are emitted and where they undergo chemical and physical processes. The air inside the University of Ghana http://ugspace.ug.edu.gh 12 box is assumed to have a homogenous concentration. The model uses that assumption to estimate the average pollutant concentration anywhere within the air shed. The following equation represents the Box model [7]: 𝑑𝐶 𝑑𝑡 𝑣 = 𝑄𝐴 + 𝑢𝐶𝑖𝑛𝑊𝐻 − 𝐶𝑊𝐻 (2.1) where 𝑄 is the pollutant emission rate per unit area, 𝐶 is the homogenous species concentration within the air shed, 𝑣is the volume described by the box, 𝐶𝑖𝑛 represents the species concentration entering the air shed and A is the horizontal area of the box. Also, 𝑢 is the wind speed normal to the box, W is the width of the box and 𝐻, the mixing height. Although the Box Dispersion model is useful, the model is unsuitable for the modeling of particle concentrations since the model simulates the formation of pollutants within the box without providing any information on the local concentrations of the pollutants [7]. 2.1.2. Lagrangian Dispersion Model Lagrangian models are similar to Box models in a sense that the models define an air shed as a box containing an initial concentration of pollutants. However, the Lagrangian model follows the trajectory of the box as the particles move downwind. The Lagrangian model then calculates the air pollution dispersion by computing the statistics of the trajectories of a large number of the pollution plume parcels. The Lagrangian equation has the following form [7]: 𝑐(𝑟, 𝑡) = 𝑝(𝑟, 𝑡 𝑟,′ 𝑡′ 𝑆 𝑟,′ 𝑡′ 𝑑𝑟′𝑑𝑡′ 𝑡 −∞ (2.2) University of Ghana http://ugspace.ug.edu.gh 13 where 𝑐(𝑟, 𝑡) is the average pollutant concentration at location, ‗𝑟‘ at time, ‗t‘; S 𝑟,′ 𝑡′ is the source emission term and 𝑝(𝑟, 𝑡 𝑟,′ 𝑡′ , the probability function that an air parcel is moving from one location at time, 𝑡 to another location at time, 𝑡′ . The disadvantage of the Lagrangian models is that, the models are limited when results from the model‘s prediction are compared with actual measurements, because measurements are made at stationary points, while the model predicts pollutants concentration based upon a moving reference grid [8]. 2.1.3. Eulerian Dispersion Model The Eulerian model is similar to the Lagrangian model but the model tracks the movement of a large number of pollution plume parcels as the particles move from their initial location. However, the models differ in a sense that the Eulerian model uses a fixed three dimensional Cartesian grid as a frame of reference rather than a moving frame of reference. The Eulerian models solve an equation of conservation of mass for a given pollutant. The equation generally follows the following form [8]: 𝜕 𝑐𝑖 𝜕𝑥 = −𝑈 .∇ 𝑐𝑖 − ∇U. 𝑐𝑖𝑈 ′ + 𝐷∇ 2 𝑐𝑖 + 𝑆𝑖 (2.3) where 𝑈 the wind field vector; 𝑈 the average wind field vector; 𝑈 ′ is the fluctuating wind fields vector; 𝑐𝑖 is the average pollutant concentration; 𝐷 the molecular diffusivity and Si is the source term. The term −𝑈 .∇ 𝑐𝑖 is hyperbolic, the turbulent diffusion is parabolic and the source term is generally defined by a set of differential equations making it difficult to solve. University of Ghana http://ugspace.ug.edu.gh 14 2.1.4. Gaussian Dispersion Model Gaussian type models are the most common dispersion models used in atmospheric dispersion modeling. The term ―Gaussian‖ refers to the statistical concept in which a group of arranged values follows a bell-shaped curve distribution. This type of model assumes that the pollutant disperses according to the normal statistical distribution. At the point of release, the pollutant concentration is at maximum and decreases in both lateral and vertical directions following the normal distribution. The analytical and numerical models used in this comparative study were developed based on Gaussian equation which is used for point source emissions in general [8]: C x, y, z, H = Q 2πuςyςz exp − y2 2ςy 2 exp z−H 2 2ςz 2 exp z+H 2 2ςz 2 (2.4) where 𝐶 is the steady-state concentration at point (𝑥,𝑦, 𝑧) in µg/m3; 𝑄 is pollutant emission rate;𝜍𝑦and 𝜍𝑧 are the standard deviations of lateral and vertical spread parameters; 𝑦 the horizontal distance from plume centerline; 𝑧 the vertical distance from ground level and 𝐻 is the effective stack height. The first exponential term represents the lateral dispersion and vertical dispersion is described by the second exponential term. The terms 𝜍𝑦 and 𝜍𝑧 in the equation represent the standard deviation of the horizontal and vertical distributions of the plume of the pollutant [8]. High standard deviation values would result from an unstable, turbulent atmosphere, whereas low values would occur in less turbulent atmospheric conditions. In older models, these coefficients are defined by stability classes created by Pasqual and Gifford and they increase as the downwind distance increases [8]. University of Ghana http://ugspace.ug.edu.gh 15 The Gaussian model is based on the following assumption:  The emission must be constant and uniform; the wind direction and speed are constant;  Net downwind diffusion is negligible compared to vertical and crosswind diffusion;  The terrain is relatively flat; there is no deposition or absorption of the pollutant and the vertical and crosswind diffusion of the pollutant follow a Gaussian distribution [8]. 2.2. SYSTEM OF RADIOLOGICAL PROTECTION Radionuclide substances are produced in small quantities but are potentially so detrimental that the presence of the substances in the environment even at low levels is sufficient to cause concern. Some of which are natural components of the environment but many of which are present largely as a result of human activity [9]. Sources of anthropogenic radionuclide substance can be related to nuclear accidents. Despite many safety features propagated by supervisory bodies in the nuclear field, an accident did occur in 1979 at the Three Mile Island PWR near Harrisburg, Pennsylvania. A maintenance error and a defective valve led to a loss-of-coolant accident [9, 10]. On March 11, 2011 a strong earthquake followed by high tsunami and fires damaged three reactors and a fuel pond at the Fukushima Dai-ichi Nuclear Power Plant in Japan with releases of radionuclide substances to the atmosphere and the sea. According to the NISA report, from 1.3 × 107𝐵𝑞 to 1.5 × 107𝐵𝑞 of 131I and about 1.3 × 1015𝐵𝑞 to 1.6 × 1016𝐵𝑞 of 137C were released to the atmosphere and the sea. The consequences of this accident at the beginning estimated as level 4 were raised to University of Ghana http://ugspace.ug.edu.gh 16 the maximum level of 7 on the INES, although the amount of discharged radioactive materials comprised approximately 10% of the Chernobyl accident only [11]. A system of radiological protection has been defined by the International Commission on Radiological Protection (ICRP), an international organization that brings out various publications on radiation protection. The system of radiological protection involves three fundamental principles; justification, optimization and dose limitation, [12]. 2.3. NUMERICAL PLUME MODELS According to McRae and Gregory [13], the plume is broken up into independent elements (plume segments or sections) whose initial features and time dynamics are a function of emission conditions and local meteorological conditions encountered by the plume segment. This model treats time varying transport conditions and changes in wind direction and speeds. Segments are sections of a Gaussian plume. Therefore, the concentration of air pollution for each segment can be calculated using the Gaussian plume model. The trajectory plume model then adds segments together to get concentration distributions over the entire plume. Arkhipov et al. [14] studied numerical modeling of pollutant dispersion by the stochastic discrete particles method and concluded that the comparison of numerical modeling and analytical results demonstrates a good fit. The stochastic discrete particles method for turbulent pollutant dispersion is often used in practice. They proposed a special procedure of random discrete particles walks which results in a solution that satisfies Richardson's law. Both approaches are found to have different solutions, thus the use of variances of random velocity for the construction of dispersion of random walking dispersion is University of Ghana http://ugspace.ug.edu.gh 17 not adequate for the turbulent diffusion process. They have shown that the suggested procedure gives a good solution that is coincident with a self-similarity solution and a one-dimensional finite-difference solution in both instantaneous and continuous discharges. Noel [15] put forward that little seems to have been published which gives useful guidelines concerning the accuracy of approximations used in air pollution transport models. In fact, he says ―In general, the application of plume models requires numerical solutions of the governing conservation equations which are highly non- linear. Rather than obtain such ‗exact‘ solutions it is easier to use various approximations to them. Sometimes, the accuracy of such approximations relative to the exact solution is quite uncertain. Noel‘s research calculates numerical solutions using the plume rise model of Hoult, Fay, and Forney [15]. Various asymptotic approximations are assessed over a wide range of parameters. Mass, momentum, energy, turbulence, temperature, and wind speed and direction are the parameters used in these approximations. Osalu et al. [16] states that ―turbulence is a meteorological quantity than can only be approximated and then only in the most ideal circumstances‖. Meteorological conditions are inherently variable and this variability causes an uncertainty in the accuracy of meteorological parameters. Plume models which incorporate these parameters inherit this uncertainty. Therefore, numerical and approximate solutions of air pollution models are dependent on the limits of meteorological modeling used to calculate the parameters for pollution models. University of Ghana http://ugspace.ug.edu.gh 18 2.4. ADVECTION-DIFFUSION EQUATION SOLUTION METHODS Many methods are described in literature which numerically solves partial differential equations of ADE. In particular, Reeden [17] compared several algorithms and variations of each which are used for solving the two-dimensional Advection Equation: ∂c ∂t + u ∂c ∂x + v ∂c ∂y = 0 (2.5) where 𝑐 is the concentration, and 𝑢 and 𝑣 are the 𝑥 𝑎𝑛𝑑 𝑦 components of the wind speed. 2.4.1. Diffusion Equation Solution Methods There are many technical techniques for solving the atmospheric diffusion equation, [24]: ∂y ∂t = ∂2c ∂x2 + ky ∂2c ∂y2 + kz ∂2c ∂z2 (2.6) Where 𝑐 is the concentration of the pollutant in question and 𝑘𝑦 𝑎𝑛𝑑 𝑘𝑧 are the 𝑦 𝑎𝑛𝑑 𝑧 diffusivity respectively. Venkatram, Goodin and Seinfeld [18] examined techniques for solving both the Advection Equation and the Diffusion Equation in their one-dimensional form in ―Numerical Solution of the Atmospheric Diffusion Equation for Chemically Reacting Flows‖. McKibbin, et al [19] studied the dispersion of heavier-than-air releases and they found that solid and liquid particulate materials that are ejected by hydrothermal or volcanic eruptions are subsequently dispersed by atmospheric wind currents while falling under gravity to the Earth‘s surface. They further noted that particle sizes are not University of Ghana http://ugspace.ug.edu.gh 19 uniform and may indeed change during flight owing to coalescence or fragmentation. Wind conditions (speed, direction and turbulence) may also change with elevation (and with time). In their attempt to determine the most important physical factors, they outlined a quantitative model that reflects the above influences on particle dispersal. Schonar [20] used a two dimensional convection diffusion equation to model air pollution. He found that the pollution dispersion is influenced by the model parameters namely diffusion coefficient, drift motion and reaction coefficient. Schonar concluded that the eventual chemical reaction between the pollutant and one or more air components acts as a pollution sink or fall out and hampers the expansion of polluted zone. Indra et al. [21], studied the classification of air pollution dispersion models and found that dispersion and emission of pollutants into the air is controlled by the prevailing meteorological conditions like wind profile, temperature profile and stability of the atmosphere. Air pollution dispersion models attempt to express the interrelationships of these factors in terms of mathematical equations. They further showed that air pollution models are based on the theories of atmospheric physics and thermodynamics. This inquest looks into the theory of air pollution modeling and critically appraises the important air pollution dispersion models. Indra and Pascal explained major types of air pollution models along with their advantages and disadvantages and concluded that, to determine the exact wind-speed and atmospheric stability class at the plume centerline height requires: (a) The prediction of the exact plume rise; (b) The exact relation between wind-speed and altitude. University of Ghana http://ugspace.ug.edu.gh 20 Indra and Pascal noted that the Gaussian models assume an ideal steady-state of constant meteorological conditions over long distances, idealized plume geometry, uniform at terrain, complete conservation of mass, and exact Gaussian distribution. Perry [22] developed a two dimensional atmospheric dispersion model for computation of the ambient air concentration of reactive pollutants emitted from ground level sources. The result showed that atmospheric chemical reactions are the most complicated and stiff part of pollutants dispersion equations. Coupling them with other physical transport processes to assemble an integrated dispersion model is a time consuming and complicated matter. Mechanism of reactions can be present in form of ordinary differential equations. Different pollutants, may present different variability characteristic due to their specific emission patterns, rates of diffusion, and transport and transformation behaviors like atmospheric reactions. Burden and Douglas [23] discussed the desirable features of the three-dimensional numerical plume models, and provide a good introduction to future work that may be undertaken. Kanevce [24] modeling has concentrated on a second-order turbulence closure for a three-dimensional boundary layer simulation with a passive pollutant. He discovered that the depth of a pollutant sink depends on the wind speed. Draxler‘s [25]modeling begins with a second-order closure to the passive pollutant transport equation and the turbulent diffusion is modeled with a zero equation model, although many other important features such as rainout, wet and dry deposition, and surface terrain have been added. University of Ghana http://ugspace.ug.edu.gh 21 2.5. AIR DISPERSION MODELING OF POLLUTANTS IN THE ATMOSPHERE Dispersion modeling uses mathematical equations describing the atmosphere, dispersion, chemical and physical processes influencing a pollutant released from sources of a given geometry to calculate concentrations at various receptors as a result of the release [26]. Dispersion models require two types of data inputs: information on the source or sources including pollutant emission rates, and meteorological data. In addition, they also need information on the topography of the study area, Figure 2.1. The models then use this information to simulate mathematically the pollutant‘s transport and dispersion. The output is air pollutant concentrations, for a particular time period, usually at specific receptor locations. A simple example of pollutants dispersion in the atmosphere is through molecular diffusion, when matter moves from a region of high concentration to a region of low concentration. However, apart from molecular diffusion, plumes spread due to other complex processes. These processes are mechanically and thermally generated turbulence and wind fluctuations [26]. Figure 2.1: A schematic representation of the input-output of an air dispersion model [26]. Source Emissions Rate and Geometry Dispersion Model Pollutant concentratio n Meteorology Topography University of Ghana http://ugspace.ug.edu.gh 22 2.5.1. Turbulence Molecules of pollutants in the air are transported from one point to another by means of turbulence. Turbulence is defined as a collective random motion involving a group of many molecules. Turbulence is made up of both thermal and mechanical eddies. Eddies are macroscopic random fluctuations from the ―average‖ flow. These turbulent eddies are responsible for the dispersion of pollutants in the atmosphere. Eddies disperse pollutants by intercepting the plume, replacing a batch of concentrated pollutants in a plume with a batch of clean air from a distance away from the plume, consequently diluting the plume and spreading it in both vertical and lateral directions [27]. 2.5.2. Mechanical Turbulence Mechanical turbulence is created through the interaction between the horizontal force exerted by one layer on an adjacent layer and the gradient of the mean velocity with height. The stronger the wind or the larger the roughness elements, the greater the mechanical turbulence hence rough surfaces such as forests or trees produce more eddies than smooth surface such as ice. Buildings, trees and other obstacles increase mechanical turbulence because these obstacles increase the horizontal forces that slow down the mean wind [28]. 2.5.3. Thermal Turbulence The thermal energy generated from the sun is absorbed by the ground. The absorbed heat is transferred into the lower atmosphere by means of conduction and/or convection thus generating thermal eddies. More eddies are created when there is strong insulation than when the energy from the sun is weak [29]. University of Ghana http://ugspace.ug.edu.gh 23 2.6. EFFECTS OF WIND FLUCTUATIONS ON DISPERSION Plume dispersion can also be caused by random shift in the wind. Pollutant concentrations are measured over a certain period of time called averaging time, for example, an averaging time of an hour. The wind direction and speed change during this period and more or less pollutant is blown towards the receptor. As a result, these random fluctuations cause the spread of the plume over a large area downwind of the source. As the plume travels downwind of the source, the pollutant spreads further in the y and z directions and the maximum concentration eventually decreases [30]. 2.7. FACTORS AFFECTING ATMOSPHERIC DISPERSION OF POLLUTANTS The factors that can affect the dispersion of pollutants in the atmosphere include: Atmospheric Stability, Meteorology, Topography and Surface Roughness [31]. 2.7.1. Atmospheric Stability Atmospheric stability is defined as the atmospheric tendency to resist or enhance vertical motion or alternatively suppress or augment existing turbulence [31]. Atmospheric stability influences the vertical movement of particles in the atmosphere which is also influenced by the temperature effect of the air. Pasquill [31] introduced a method of estimating atmospheric stability accounting for both mechanical and thermal turbulence. Atmospheric stability was classified into six categories ranging from A (very unstable) to F (very stable), Table 2.1.The categories were developed based on the wind speed, solar radiation (daytime) and cloud cover (at night). Strong insulation leads to the heating of the ground increasing the temperature of the lower part of the atmosphere, creating an unstable condition. If the wind speed rises, the University of Ghana http://ugspace.ug.edu.gh 24 vertical mechanical mixing becomes stronger than the buoyancy effects, leading to neutral stability. During the night the ground cools creating stable conditions [31]. Table 2.1: Atmospheric Stability Classification Developed by Pasquill [31] Day Night Wind speed at 10m, incoming solar Thin m/s radiation cloud <3/8 Strong Moderate slight or >3/8 cloud >2 A A-B B 2-3 A-B B C E F 3-5 B B-C C D E 5-6 C C-D D DD >6 C D DDD Overcast conditions yield D both day and night 2.7.2. Weather Conditions Meteorology is a vital element of atmospheric dispersion modeling because it determines the diluting effects of the atmosphere. The dispersion, transformation and removal of pollutants in the atmosphere depend on the meteorological conditions of the site. Hence, good and appropriate meteorological data preferably from a weather station within the area of interest is needed in order to achieve the best results from modeling. The important meteorological data needed for modeling are: temperature, wind speed, wind direction, cloud cover and atmospheric stability [32]. University of Ghana http://ugspace.ug.edu.gh 25 2.7.2.1. Air Temperature Temperature affects the buoyancy of the plume since the higher the temperature difference between ambient air and the plume, the higher the plume will rise. This in turn reduces the ground level impact of pollutants. Temperature is also important for the development of the mixing and inversion layer [32]. 2.7.2.2. Wind Speed Wind speed is one of the important meteorological parameters in dispersion modeling. Wind speed influences initial dilution of the plume leaving a source, hence the stronger the wind speed, the more rapid the dilution of the pollutants and thus the lower the concentrations at the ground level and vice versa. Mechanical turbulence that increases mixing and dilution is created by the wind and the higher the wind speed the stronger the mechanical turbulence [32]. 2.7.2.3. Wind Direction The wind rose shown in Figure 2.2 determines the predominant wind which will affect the transport pollutants in the atmosphere. The wind direction played a major role since the wind direction was used to determine the monitoring point where the monitor was set. Only monitoring points or receptors in the North-Eastern direction of the source are affected by the plume emitted. Wind direction together with other meteorological parameters determines the spatial pattern of average ground level concentration [32]. University of Ghana http://ugspace.ug.edu.gh 26 Figure 2.2: Wind Rose for Accra 2.7.3. Topography Topography also influences the dispersion of air pollutants. The term ―topography‖ refers to the surface features of land, including the configuration and elevation of both man-made and natural feature. Topographical features may impede the dispersion of pollutants, especially when the pollutants are released in low-lying areas. Surface roughness, buildings, hills, trees and obstructions are some of the topographical features that can affect pollutant dispersion in the atmosphere. The effect of surface roughness on dispersion is further discussed briefly in the next paragraph [33]. University of Ghana http://ugspace.ug.edu.gh 27 2.8. DISPERSION PARAMETERS The horizontal dispersion parameter (𝜍𝑦) used in the dispersion model was based on the results of a large number of dispersion experiments. The experiments, which were conducted over relatively flat terrain, typically involved tracer releases ranging from about 1 𝑚/𝑠 to 7.5 𝑚/𝑠 with ground-level concentration measurements at distances ranging from 100m to several kilometers. Only a few direct measurements of vertical dispersion parameters (𝜍𝑧) were made. Consequently, vertical dispersion parameters have been estimated with dispersion models using measured values of the horizontal dispersion parameter and measured concentrations. Dispersion parameter (Wind speed, Diffusion, Concentration gradient or density gradient, Emission rate, Settling Velocity, Viscosity of air) values used in the Gaussian Plume model have been summarized in Table 2.2, [34]. The basic dispersion parameter relationships are [34]: σy x = ayx 0.9031 (2.7) and ςy x = az x . x bz (x ) + cz(x) (2.8) where 𝑥 the distance from the released point, in meters; 𝑎𝑦 is a function of stability class and 𝑎𝑧 ,𝑏𝑧 ,𝑎𝑛𝑑 𝑐𝑧 are the functions of stability class and distance. The values 𝑎𝑧 ,𝑏𝑧 ,𝑐𝑧 and 0.9031 are empirical constant evaluated by fitting curves. Of these constants, 0.9031 and 𝑏𝑧are dimensionless; 𝑐𝑧 has dimensions of meters; and 𝑎𝑦 and 𝑎𝑧 have dimensions of𝑚 1−0.9031and𝑚1−𝑏 , respectively. University of Ghana http://ugspace.ug.edu.gh 28 Table 2.2: Parameter Values [Cs-137] Used in the Gaussian Plume Model [34] Parameter Symbol Value Units 𝜍𝑦 18.6 𝜍𝑧 13.1 Wind speed 𝑢 2.5-10 m/s Source heights 𝐻𝑠 [0,5] m Emission rate 𝑄𝑠 10 mg/s Diffusion parameter a 0.44 𝑚2−𝑏 Diffusion parameter b 0.76 𝑚2−𝑏 Settling velocity 𝑤𝑠𝑒𝑡 2.7e-3 m/s Deposition velocity 𝑤𝑑𝑒𝑝 5.0e-3 m/s Viscosity of air 𝜇 1.8e-5 mg/ms Gravitational 𝑔 9.8 m/s2 Particle density 𝜌 432 mg/m3 Particle radius 𝑅 1.3e-5 m 2.9. FINITE DIFFERENCE METHOD FOR SOLVING THE ADVECTION- DIFFUSION EQUATION (ADE) The finite difference method is one of several techniques for obtaining numerical solutions to the Advection Diffusion Equation. In all numerical solutions, the continuous partial differential equation (PDE) is replaced with a discrete approximation. In this context the word ‗discrete‘ means that the numerical solution is known only at a finite number of points in the physical domain. The number of those points can be selected by the user of the numerical method. In general, increasing the number of points not only increases the resolution, but also the accuracy of the numerical solution [35]. The mesh is the set of locations where the discrete solution is computed. These points are called nodes, and if one were to draw lines between adjacent nodes in the domain University of Ghana http://ugspace.ug.edu.gh 29 the resulting image would resemble a net or mesh. Two key parameters of the mesh are ∆x, ∆y, the local distances between adjacent points in space with their respective coordinates, and ∆t, the local distance between adjacent time steps. The core idea of the finite difference method is to replace continuous derivatives with so called difference formulas that involve only the discrete values associated with positions on the mesh [35]. Applying the finite difference method to a differential equation involves replacing all derivatives with difference formulas. In the advection diffusion equation there are derivatives with respect to time, and derivatives with respect to space. By using different combinations of mesh points in the difference formulas results in different schemes. In the limit as the mesh spacing (∆x, ∆y and ∆t) go to zero, the numerical solution obtained with any useful1 scheme will approach the true solution to the original differential equation. However, the rate at which the numerical solution approaches the true solution varies with the scheme. In addition, there are some practically useful schemes that can fail to yield a solution for bad combinations of ∆x, ∆y and ∆t [35]. The finite difference method involves using discrete approximations like: ∂C ∂x ≈ Ci+1 − Ci ∆x where the quantities on the right hand side are defined on the finite difference mesh. Approximations to the governing differential equation are obtained by replacing all continuous derivatives by discrete formulas [35]. University of Ghana http://ugspace.ug.edu.gh 30 CHAPTER THREE RESEARCH METHODOLOGY The research examined the Gaussian Plume model and numerical schemes adopted to solve the steady state Advection-Diffusion Equation. Air pollution transport models and the partial differential equations used in the air pollution modeling were investigated. The analytical and numerical solutions of the mathematical equation of the Advection-Diffusion Equation were solved by using codes developed in Matrix Laboratory (MATLAB) and Microsoft Office Excel spread sheet. The MATLAB codes showed in Appendices II—V were implemented and the results obtained were presented and discussed in Chapter Four. From the code implementation, the concentration of Cs-137 radionuclide at a given location, the deposition and mass of the radionuclide deposited at each receptor point would be determined. 3.1. PHYSICAL DESCRIPTION OF THE TRANSPORT PHENOMENON When a plume is discharged from an elevated source, the plume would spread vertically until the plume‘s lower edge reaches the ground. Figure 3.1 shows the Gaussian-plume approach to the atmospheric dispersion modeling which describes the bell-shaped (Gaussian) distribution of concentrations in the horizontal and vertical directions. The Figure also illustrates the plume rise and the normal (Gaussian) distribution of the contaminant concentrations in the horizontal and vertical directions. As the standard deviation of the vertical distribution (𝜍𝑧) increases, the plume will make a significant contribution of Cs-137 concentrations on the ground [36]. University of Ghana http://ugspace.ug.edu.gh 31 Similarly, if strong inversion layer is located at some height (𝑧𝑖) above the source, then the plume will have difficulty expanding vertically and effectively be trapped between the inversion and the ground, [36]. The Gaussian-plume model is derived assuming ‗steady-state‘ conditions, i.e. the Gaussian-plume dispersion model does not depend on time, although the model does represent an ensemble of time average. The meteorological conditions are assumed to remain constant during the dispersion from source to receptor, which is effectively instantaneous. Figure 3.1: Gaussian Shape of the Plume Cross-Sections Relative to the Plume Centerline [38]. University of Ghana http://ugspace.ug.edu.gh 32 3.2. DERIVATION OF THE ADVECTION-DIFFUSION EQUATION In order to derive an equation describing the distribution of mass of Cs-137 within the plume, the transport of mass within a small control volume was considered as shown in Figure 3.2. Figure 3.2: Dispersion Control Volume of the ADE University of Ghana http://ugspace.ug.edu.gh 33 The Gaussian plume equation is a particular solution to the ADE under the following assumptions: (i) State conditions, 𝜕𝐶 𝜕𝑡 = 0 (ii) Steady Constant wind speed at height,𝑕 (𝑢 does not depend on 𝑧) (iii)Constant eddy diffusivity (𝐾 does not depend on 𝑦 𝑎𝑛𝑑 𝑧) (iv) Mass is conserved , 𝐶𝑑𝑦𝑑𝑧 = 𝑄 ; 𝑓𝑜𝑟 𝑥 > 0 ∞ −∞ ∞ −∞ . By considering Figure 3.2; Net rate of change of mass flow = (Mass flow rate in) — (Mass flow rate out). Hence, 1. Changes in 𝑥 − 𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛 due to ―Advection‖ (i.e. Transport by the wind, u and mean velocity,𝑣). Mass flow rate in = 𝐶𝑢𝐴𝑦𝑧 𝜇𝑔/𝑚 3 𝑚/𝑠 𝑚2 (3.1) Mass flow rate out = 𝐶𝑢𝐴𝑦𝑧 + 𝜕 𝜕𝑥 𝐶𝑢𝐴𝑦𝑧 𝑑𝑥 (3.2) Net Rate of Change= − 𝜕 𝜕𝑥 𝐶𝑢𝐴𝑦𝑧 𝑑𝑥 = − 𝜕 𝜕𝑥 𝐶𝑢 𝑉 (3.3) Since 𝑉 = 𝑑𝑥𝑑𝑦𝑑𝑧; 𝐴𝑥𝑦 = 𝑑𝑥𝑑𝑦 ;𝐴𝑦𝑧 = 𝑑𝑦𝑑𝑧 and 𝐴𝑥𝑧 = 𝑑𝑥𝑑𝑧. 2. Changes in 𝑦 − 𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛 via ―Turbulent Diffusion‖ Mass flow rate in = −𝐴𝑥𝑧 𝜕 𝜕𝑦 𝐾𝑦𝐶 𝑑𝑦 𝑚2 𝑚−1 𝜇𝑔𝑚−3 𝑚2/𝑠𝑒𝑐−1 (3.4) Net rate of change = ∂ ∂y ∂ ∂y KyC V (3.5) 3. Changes in 𝑧 − 𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛 via ―Turbulent Diffusion‖ University of Ghana http://ugspace.ug.edu.gh 34 Mass flow in = −𝐴𝑥𝑦 ∂ ∂z KzC dz m2 m−1 μgm−3 m2/sec−1 (3.6) Net rate of change = ∂ ∂z ∂ ∂z KyC V (3.7) where 𝑅𝑎𝑡𝑒 = 𝐾 𝛿𝑐 𝛿𝑥 and K is termed an ―eddy diffusivity‖ with units of 𝑚2/𝑠𝑒𝑐. 4. Given that the net rate of change in Control Volume 𝑉 𝛿𝑐 𝛿𝑡 is the change in all three directions, an overall expression in terms of 𝑥, 𝑦 𝑎𝑛𝑑 𝑧 as: [i.e. the sum of equations (3.3), (3.5) and (3.7)], and the source term is the partial differential along the 𝑥, 𝑦 𝑎𝑛𝑑 𝑧 components. ∂C ∂t = − U ∂C ∂x + ∂ ∂y Ky ∂C ∂y + ∂ ∂z Kz ∂C ∂z + Qδ x δ y δ(z− H) (3.8) where 𝜕𝐶 𝜕𝑡 = 0 is the steady-state; 𝑈𝜕𝐶 𝜕𝑥 is the advection (i.e. Transport by the mean wind, u); 𝐾 is the eddy diffusivity representing the intensity of turbulent motions and varying with stability; where 𝜕 𝜕𝑦 𝐾𝑦 𝜕𝐶 𝜕𝑦 + 𝜕 𝜕𝑧 𝐾𝑧 𝜕𝐶 𝜕𝑧 is the effect of turbulent diffusivity and 𝑆 is the source term. The equation (3.8) is known as the ―Advection Diffusion Equation‖. 3.3. ANALYTICAL SOLUTION OF THE ADVECTION-DIFFUSION EQUATION 3.3.1. Assumptions and Boundary Conditions [38]; A number of simplifying assumptions were made to permit the derivation of the analytical solution of the advection-diffusion equation. A1. The contaminant was emitted at a constant rate Q [mg/s] from a single point source𝑥 = 0,0,𝐻 located at height H above the ground surface. The source term was written as; 𝑆 𝑥 = 𝑄𝛿 𝑥 𝛿 𝑦 𝛿 𝑧 − 𝐻 [37,38]. University of Ghana http://ugspace.ug.edu.gh 35 where 𝛿(. ) is the Dirac delta function. The units of the delta function are in [m-1]. For the stack-like configuration from Figure 3.1, the height was actually an effective height ( 𝐻 = 𝑕 + 𝛿𝑕), which is the sum of the actual stack height 𝑕 and the plume rise 𝛿𝑕 that arose from buoyant effects. A2. Constant wind speed at height,𝑕 (u did not depend on z) and aligned with the positive 𝑥-axis so that 𝑢 = (𝑢, 0,0) for some constant 𝑢 ≥ 0. A3. The solution was for steady state, which was reasonable if the wind velocity and all other parameters were independent of time and the time scale of interest was long enough, where at steady-state conditions, 𝜕𝐶 𝜕𝑡 = 0. A4. Constant eddy diffusivity (𝐾 did not depend on y or z). The eddy diffusivities were functions of the downwind distance,𝒙 only, and diffusion was isotropic so that 𝐾𝑥 𝑥 = 𝐾𝑦 𝑥 = 𝐾𝑧 𝑥 =:𝐾(𝑥), defined by: ςy 2 2Ky x u and ςz 2 2Kz x u . A5. The wind velocity was sufficiently large that diffusion in the 𝑥-direction was much smaller than advection; then the term 𝐾𝑥𝛿𝑥 2𝐶 was neglected, i.e. 𝑘 𝑢 ≈ 1. A6. Variations in topography were negligible so that the ground surface was taken as the plane, 𝑧 = 0. A7. The contaminant did not penetrate the ground. A8. Mass was conserved 𝐶𝑑𝑦𝑑𝑧 = 𝑄 𝑓𝑜𝑟 𝑥 > 0 ∞ −∞ ∞ −∞ 𝜇𝑔 . University of Ghana http://ugspace.ug.edu.gh 36 3.3.2. REDUCTION OF ADE TO BOUNDARY VALUE PROBLEM (BVP) TRANSFORMATION By applying A3 (steady state) to Equation (3.8), the Advection Diffusion Equation was reduced to: u ∂C ∂x = K ∂2C ∂y2 + ∂2C ∂z2 + Qδ x δ y δ(z− H) (3.9) Integrating both sides of equation (3.9) over the integral 𝑥 ∈ −𝑥, 𝑥 ; 𝑢 𝜕𝐶 𝜕𝑥 𝑑 −𝑑 = 𝐾 𝜕2𝐶 𝜕𝑦2 𝑑 −𝑑 + 𝜕2𝐶 𝜕𝑧2 + 𝑄𝛿 𝑥 𝛿 𝑦 𝛿(𝑧 − 𝐻) 𝑑 −𝑑 𝑑 −𝑑 (3.10) 𝑢 𝐶 𝑑,𝑦, 𝑧 − 𝐶 −𝑑,𝑦, 𝑧 = 2𝑑𝐾 𝜕2𝐶 𝑑 𝜕𝑦2 + 𝜕2𝐶 𝑑 𝜕𝑧2 + 𝑄𝛿 𝑥 𝛿 𝑦 𝛿(𝑧 − 𝐻 𝑑 −𝑑 (3.11) where 𝐶 𝑑 𝑦, 𝑧 = 1 2𝑑 𝐶 𝑥,𝑦, 𝑧 𝑑𝑥 𝑑 −𝑑 represents the average value of concentration over the interval – 𝑑,𝑑 and 𝛿 𝑥 𝑑𝑥 𝑑 −𝑑 = 1. Now, limit as 𝑑 → 0+, equation (3.11) reduced to: 𝐶 0,𝑦, 𝑧 = 𝑄 𝑢 𝛿 𝑦 𝛿(𝑧 − 𝐻) (by making use of the assumption that 𝐶 −𝑑,𝑦, 𝑧 = 0). Solving the PDE in equation (3.9) on the interval 𝑥 > 0, the delta function term can be neglected and the PDE reduced to BVP ( Boundary Value Problem) given by: u ∂C ∂x = K ∂2C ∂y2 + ∂2C ∂z2 (3.12) with boundary conditions given as: University of Ghana http://ugspace.ug.edu.gh 37 C 0, y, z = Q u δ y δ z− H ; 𝐶 ∞,𝑦, 𝑧 = 0; 𝐶 𝑥, ±∞, 𝑧 = 0; 𝐶(𝑥,𝑦,∞) = 0; 𝐾 𝜕𝐶 𝜕𝑧 𝑥,𝑦, 𝑜 = 0. Replacing 𝑥 in equation (3.12) with new independent variable (change of variables) 𝑟 = 1 𝑢 𝐾 𝑥 0 𝑥ı 𝑑𝑥ı under the derivatives transform, according to 𝜕𝑥 = 𝐾 𝑢 𝜕𝑟 and defining 𝑐 𝑟,𝑦, 𝑧 = 𝐶(𝑥, 𝑦, 𝑧), the transformed problem is given by [38]: 𝜕𝐶 𝜕𝑟 = 𝜕2𝐶 𝜕𝑦2 + 𝜕2𝐶 𝜕𝑧2 (3.13) with boundary conditions: 𝑐 0,𝑦, 𝑧 = 𝑄 𝑢 𝛿 𝑦 𝛿 𝑧 − 𝐻 ; 𝑐(𝑟,𝑦,∞) = 0 𝑐 ∞,𝑦, 𝑧 = 0; and 𝜕𝐶 𝜕𝑧 𝑟, 𝑦, 𝑜 = 0 Assuming a separable solution of the form (i.e. the dependence of the solution on y and z can be separated according to): 𝑐 𝑥,𝑦, 𝑧 = 𝑄 𝑢 𝑎 𝑥,𝑦 .𝑏(𝑥, 𝑧) (3.14) By substituting the separable solution into the PDE in Equation (3.13) to obtain: 𝑎 𝜕𝑏 𝜕𝑥 + 𝑏 𝜕𝑎 𝜕𝑥 = 𝑏 𝜕2𝑎 𝜕𝑦2 + 𝑎 𝜕2𝑏 𝜕𝑧2 (3.15) Equation (3.15) is reduced to two-dimensional diffusion equation as: 𝜕𝑎 𝜕𝑥 = 𝜕2𝑎 𝜕𝑦2 (3.16) 𝑓𝑜𝑟 0 ≤ 𝑥 < ∞ 𝑎𝑛𝑑 −∞ < 𝑦 < ∞ with boundary conditions: University of Ghana http://ugspace.ug.edu.gh 38 𝑎 0,𝑦 = 𝛿 𝑦 , 𝑎 ∞,𝑦 = 0, 𝑎𝑛𝑑 𝑎 𝑥, ±∞ = 0 𝜕𝑏 𝜕𝑥 = 𝜕2𝑏 𝜕𝑧2 (3.17) 𝑓𝑜𝑟 0 ≤ 𝑥 < ∞ 𝑎𝑛𝑑 −∞ < 𝑧 < ∞, with boundary conditions: 𝑏 0, 𝑧 = 𝛿 𝑧 − 𝐻 , 𝑏 ∞, 𝑧 = 0, 𝑎𝑛𝑑 𝜕𝑏 𝜕𝑧 𝑥, 0 = 0 In both problems, the variable ′𝑥′ can be viewed as a time-like variable and so the boundary at𝑥 = 0 (which contains the delta function) acts as initial conditions for the respective diffusion equations. 3.3.3. SOLUTION OFTRANSFORMED ADVECTION-DIFFUSION EQUATION By considering equation (3.14) [38]: 𝑐 𝑥, 𝑦, 𝑧 = 𝑄 𝑢 𝑎 𝑥,𝑦 . 𝑏(𝑥, 𝑧) For 𝑎 𝑥,𝑦 the Laplace transform of the PDE in 𝑥 is: 𝜌𝑎 − 𝑎 0,𝑦 = 𝜕2𝑎 𝜕𝑦2 (3.18) Where 𝑎 𝜌,𝑦 = ℒ𝑥 𝑎(𝑥,𝑦) = 𝑒 −𝜌𝑥𝑎 𝑥,𝑦 𝑑𝑥 ∞ 0 and 𝜌 is the transform variable. By applying the source boundary condition, equation (3.18) reduces to: 𝜕2𝑎 𝜕𝑦2 − 𝜌𝑎 = −𝛿 𝑦 (3.19) By taking the Laplace transform in 𝑦 as [38]: 𝜂2𝑎 − 𝜂𝑎 𝜌, 0 − 𝜕𝑎 𝜕𝑦 𝜌, 0 − 𝜌𝑎 = −1 (3.20) University of Ghana http://ugspace.ug.edu.gh 39 where 𝑎 𝜌, 𝜂 = ℒ𝑦 𝑎 (𝜌, 𝑦) = 𝑒 −𝜂𝑦 𝑎 𝜌,𝑦 𝑑𝑦 ∞ 0 and 𝜂 is the transform variable. Equation (3.20) can be simplified to obtain: 𝑎 𝜌, 𝜂 = 𝜂𝐶1+𝐶2 𝜂2−𝜌 (3.21) where 𝐶1 = 𝑎 𝜌, 0 and 𝐶2 = 𝜕𝑦𝑎 𝜌, 0 − 1 Now, by applying the inverse transform in 𝜂 to obtain: 𝑎 𝜌,𝑦 = 𝐶1 cosh 𝜌𝑦 − 𝐶2 2 𝜌 𝑒 𝜌𝑦 − 𝑒− 𝜌𝑦 (3.22) As 𝑎 → 0, 𝑦 → ∞ and 𝐶1 = 𝐶2 𝜌 then Equation( 3.22) for 𝑎 reduced to 𝑎 𝜌,𝑦 = 𝐶2 𝜌 𝑒− 𝜌𝑦 (3.23) Similarly, applying the inverse transform in 𝜌 to get 𝑎(𝑥,𝑦) = 𝐶2 𝜋𝑥 𝑒𝑥𝑝 −𝑦2 2𝑥 (3.24) Now, by employing the delta function identity [38]: 𝛿 𝑦 = lim𝑥→0 𝑒𝑥𝑝 −𝑦2 2𝑥 / 2𝜋𝑥 and 𝐶2 = 1/2 , equation (3.24) gives: 𝑎 𝑥,𝑦 = 1 2σz 2 𝑒 −𝑦2/2𝑥 (3.25) Similarly, for 𝑏(𝑥, 𝑧) in equation (3.14), the Laplace transform of the PDE in 𝑥 is given by: 𝜕2𝑏 𝜕𝑧 2 − 𝜌𝑏 = 𝛿(𝑧 − 𝐻) (3.26) University of Ghana http://ugspace.ug.edu.gh 40 where 𝑏 𝜌, 𝑧 = ℒ𝑟 𝑏(𝑥, 𝑧) . Taking the Laplace transform in 𝑧 and defining: 𝑏 𝜌, 𝜁 = ℒ𝑧 𝑏 𝜌, 𝑧 Equation (3.26) reduced to: 𝜁2𝑏 − 𝜁𝑏 𝜌, 0 − 𝜕𝑏 𝜕𝑧 𝜌, 0 − 𝜌𝑏 = −𝑒−𝜁𝐻 (3.27) Appling the transformed NBC (Neumann Boundary Condition): 𝜕𝑧𝑏 𝜌, 0 = 0 𝑏 𝜌, 𝜁 = 𝜁𝑏 𝜌 ,0 −𝑒−𝜁𝐻 𝜁2−𝜌 (3.28) Now, by taking the inverse transform in 𝜁 to obtain: 𝑏 𝜌, 𝑧 = 𝑏 𝜌, 0 𝑐𝑜𝑠𝑕 𝜌𝑧 − 1 𝜌 𝑠𝑖𝑛𝑕 𝜌(𝑧 − 𝐻) (3.29) With the condition that as𝑏 → 0, 𝑧 → ∞, 𝑏 𝜌, 0 = exp − 𝜌𝐻 / 𝜌. Hence, b ρ, z = 1 2 ρ e− ρ(z−H) + e− ρ(z+H) (3.30) The inverse transform in 𝜌 is given by: b x, z = 1 2σz 2 e −(z−H)2 2x + e−(z+H) 2 2x (3.31) By substituting equations (3.25) and (3.31) into equation (3.14) yields: c x, y, z = Q 2πuσyσz exp −y2 2σz 2 −(z−H)2 2σz 2 + exp −(z+H)2 2σz 2 (3.32)Hence, the Gaussian Plume equation (3.32) becomes [38]: University of Ghana http://ugspace.ug.edu.gh 41 C x, y, z = Q 2πuσyσz exp −(z−h)2 2σz 2 + exp −(z+h)2 2σz 2 −(y)2 2σy 2 (3.33) where 𝐶 𝑥,𝑦, 𝑧 is the pollutant concentration as a function of downwind position 𝑥,𝑦, 𝑧 ; 𝑄 is the mass emission rate; 𝑢 is the wind speed in the 𝑥 −direction evaluated at effective release height;𝜍𝑦 𝑎𝑛𝑑 𝜍𝑧are Gaussian parameters (values depend upon downwind distance, x); and 𝑕 is the effective stack height. 3.4. FINITE DIFFERENCE SOLUTION OF THE ADVECTION-DIFFUSION MODELFOR RADIONUCLIDE (Cs-137) ATMOSPHERIC TRANSPORT Finite difference approximations to higher order derivatives can be obtained with the additional manipulations of the Taylor‘s series expansion about𝐶(𝑥𝑖). C1+i = Ci + h ∂C ∂x + h2 2! ∂2C ∂x2 + h3 3! ∂3C ∂x3 +⋯ (3.34) C1−i = Ci − h ∂C ∂x + h2 2! ∂2C ∂x2 − h3 3! ∂3C ∂x3 +⋯ (3.35) Adding equations (3.34) and (3.35) yields: C1+i + C1−i = 2Ci+h 2 ∂ 2C ∂x2 + 2 h4 4! ∂4C ∂x4 +⋯ (3.36) Solving for 𝜕2𝐶 𝜕𝑥 2 gives ∂2C ∂x2 = C1+i−2C i+C i−1 h2 + 2 h4 4! ∂4C ∂x4 +⋯ (3.37) where the last term on the right hand side is a temporal truncation error and is denoted by 𝜑(𝑕4). Thus, the central differential approximation for the second order derivative is obtained by: ∂2C ∂x2 = C1+i−2C i+C i−1 h2 + φ(h4) (3.38) University of Ghana http://ugspace.ug.edu.gh 42 3.4.1. Choice of Method The Crank-Nicholson method was used for the finite difference solution because it is unconditionally stable and is second-order accurate in step size as well as in space. Also, the method used the assumption that the partial differential equation is satisfied not just at the grid points but also at points in time halfway between two grid points [39]. 3.4.2. Crank- Nicholson Solution of the ADE Consider a two dimensional function 𝐶(𝑥, 𝑦, 𝑡) which is a function for both space and time, the central difference for the first order equation is obtained by modifying equation below by considering 𝐶(𝑥,𝑦, 𝑡) as two dimensional in space. Now if we consider 𝛿𝑥 = 𝛿𝑦 = 𝑕, so the central difference for 𝜕𝐶 𝜕𝑥 ; 𝜕𝐶 𝜕𝑥 = 𝐶𝑖+1,𝑗−𝐶𝑖−1,𝑗 2𝑕 + 𝜑(𝑕2) (3.39) And for 𝜕𝐶 𝜕𝑦 is: ∂C ∂y = C i ,j+1−C i ,j−1 2h + φ(h2) (3.40) Similarly, the second order difference equations are: ∂2C ∂x2 = C i+1,j−2C i ,j +Ci−1,j h2 + φ(h4) (3.41) ∂2C ∂y2 = C i ,j+1−2C i ,j +Ci ,j−1 h2 + φ(h4) (3.42) where 𝑖 𝑎𝑛𝑑 𝑗 are index in 𝑥 𝑎𝑛𝑑 𝑦 values respectively. University of Ghana http://ugspace.ug.edu.gh 43 By Crank-Nicholson method, ∂C ∂x i,j+1 2 = C i ,j+1−C ij 2(k 2) = C i ,j+1−C ij k (3.43) where 𝑘 is the number of step size in the 𝑦 direction and 𝑕 is the mesh size. Also, ∂ 2C ∂y2 i,j+1 2 = 1 2 C i−1,j−2C i ,j +Ci+1,j h2 + C i−1,j+1−2C i ,j+1+Ci+1,j+1 h2 (3.44) Similarly, ∂ 2C ∂z2 i,j+1 2 = 1 2 C i−1,j−2C i ,j +Ci+1,j h2 + C i−1,j+1−2C i ,j+1+Ci+1,j+1 h2 (3.45) ∂C ∂x = K u ∂2C ∂y2 + ∂2C ∂z2 + S (3.46) Hence, applying the finite difference on Equation (3.46), using equations (3.43), (3.44) and (3.45) yielded the finite difference equation for 2-Dimensional ADE as: Ci,j+1 − Cij k = K u 1 2 Ci−1,j − 2Ci,j + Ci+1,j h2 + Ci−1,j+1 − 2Ci,j+1 + Ci+1,j+1 h2 + 1 2 Ci−1,j − 2Ci,j + Ci+1,j h2 + Ci−1,j+1 − 2Ci,j+1 + Ci+1,j+1 h2 + S (3.47) Let 𝐾 𝑢 = 1, with no restriction on 𝐾 𝑢 and the source term 𝑆 = 0, then Equation (3.47) reduced to: C i ,j+1−Cij k = 1 2h2 2Ci−1,j − 4Ci,j + 2Ci+1,j + 2Ci−1,j+1 − 4Ci,j+1 + 2Ci+1,j+1 (3.48) = 1 h2 Ci−1,j − 2Ci,j + Ci+1,j + Ci−1,j+1 − 2Ci,j+1 + Ci+1,j+1 University of Ghana http://ugspace.ug.edu.gh 44 Ci,j+1 − Cij = k h2 Ci−1,j − 2Ci,j + Ci+1,j + Ci−1,j+1 − 2Ci,j+1 + Ci+1,j+1 (3.49) Thus, −Cij − k h2 Ci−1,j − 2Ci,j + Ci+1,j = −Ci,j+1 + k h2 Ci−1,j+1 − 2Ci,j+1 + Ci+1,j+1 (3.50) Let 𝑘 𝑕2 = 1, with no restriction on 𝑘 𝑕2 , the difference formula yields: −Cij − Ci−1,j + 2Ci,j − Ci+1,j = −Ci,j+1 + Ci−1,j+1 − 2Ci,j+1 + Ci+1,j+1 (3.51) Hence, Ci−1,j+1 − 3Ci,j+1 + Ci+1,j+1 = −Ci−1,j + Cij − Ci+1,j (3.52) University of Ghana http://ugspace.ug.edu.gh 45 3.4.3. NUMERICAL FLOW CHART FOR THE ADE OF THE GAUSSIAN MODEL The flow chart in Figure 3.3 demonstrated the input and output of the implemented values and the conditions involve in the coding of the ADE (Equation 3.8) and the Numerical Solution (Equation 3.52). Figure 3.3: Flow Chart for the ADE of the Gaussian Plume Model input end points X and Y ; parameters- velocity vector; the source height START initialize fields: grav, mu,Wset, Wdep, diam set receptor locations and coordinates in steps of CM=CN =0, if u ∈ [-d,d] If Umin< U For𝑗 = 1,… ,𝑑 − 1, Do step 6 to step 8 compute C with new clist, clist=[1e-3,2e-3,2e-4, 4e-4] if input=Cx solve for ground level concentration End Program NO YES YES NO University of Ghana http://ugspace.ug.edu.gh 46 3.4.4. STRUCTURED PROGRAM The code configuration incorporated the setting up the data and the systematic organization of the subroutines of the parameters. Steps 1—10, gave the highlights of the processes involved in the coding of the ADE (Equation 3.8) and the Numerical Equation of the ADE (Equation 3.52), which was implemented in Appendix II-VI. Step 1 — Physical parameters for the atmospheric dispersion problem were set up. Step 2 — Specification of parameters (Dispersion rate: change in plume centerline). Step 3 — Input stacks emission source data and receptor diameter (step size, t=300s). Step 4 — Locations of receptors with deposition measurements were input. Step 5 — Sigma coefficients based on stability class were determined. Step 6 — Cut-off velocity, below which concentration = 0. Step 7 — Subroutines from data file setparams run. Step 8—Setting plotting parameters with corresponding mesh points for contour plot. Step 9—The ground level concentration and the total deposition at a set of r receptor jar locations were calculated. Step10— Comparison of the integrated crosswind was implemented. University of Ghana http://ugspace.ug.edu.gh 47 3.6. SOLUTION ALGORITHM BASED ON MATLAB CONTOURING SOFTWARE The solution algorithm which was implemented in Appendices II—VI summarized the formation and configuration of functions and contour plots used in the coding of the ADE (Equation 3.52. 1. Creation of functions % Function to compute contaminant concentration in kg/m3 function 𝐶 = (𝑥, 𝑦, 𝑧,𝐻,𝑄,𝑈) % input parameters 𝑥 = 50; % receptor locations: distance along the wind direction in meters 𝑦 = 0; % receptor locations: cross-wind direction in meters 𝑧 = 0; % receptor locations: vertical height in meters 𝑄 = 0.01; % contaminant emission rate in kg/s % Determine the sigma coefficients based on stability class 𝑎𝑦 = 0.34; 𝑏𝑦 = 0.82; 𝑎𝑧 = 0.275; 𝑏𝑧 = 0.82 𝑠𝑖𝑔𝑚𝑎𝑦 = 𝑎𝑦 ∗ 𝑎𝑏𝑠 𝑥 .𝑏𝑦 .∗ (𝑥 > 0) 𝑠𝑖𝑔𝑚𝑎𝑧 = 𝑎𝑧 ∗ 𝑎𝑏𝑠 𝑥 .𝑏𝑧 .∗ (𝑥 > 0) 𝑖𝑓 𝑈 < 𝑈𝑚𝑖𝑛 𝐶 = 0 ∗ 𝑧 𝑒𝑙𝑠𝑒 University of Ghana http://ugspace.ug.edu.gh 48 𝐶0 = 2 ∗ 𝑄/(𝑝𝑖 ∗ 𝑈 ∗ 𝐻2 ∗ 𝑒𝑥𝑝9 1 ); % maximum concentration value 2. Creation of contour % Contours created for the graphics [cs,ch]= contour (xx,zz,cc, clist); % define clist Clist=[1e-5, 1e-4, 1e-3, 1e-2, 0.025, 0.05, 0.01]; Clabel (cs) 𝑥𝑙𝑎𝑏𝑒𝑙 (x(m), ylabel(m), z(m)) text (0.05*max, 0.88*zmax, ‗wind‘) hold off 3. Creation of meshgrids % Createmeshgrids to apply boundary conditions [rr,yy]= meshgrid ([1:N]*rmax/N, [-N,N]*ymax/N); [cs,ch]=contour (xx,yy,cc,clist); Grid on, hold off for i = 1 : source.n, % Sum up ground-level Iodine concentrations from each source at all mesh initialize glc = 0;dep = 0; glc = glc + ermak( xmesh-source.x(i), ymesh-source.y(i), 0.0, ... source.z(i), source.Q(i), Uwind, Wdep, Wset ); end. University of Ghana http://ugspace.ug.edu.gh 49 In chapter four, the above algorithm was implemented and approximate numerical solutions were found. Furthermore, the boundaries of the infected region as well as the dynamics of the polluted zone were analyzed. 3.6.1. PROGRAM IMPLEMENTATION The simulated results of Equation (3.33) and Equation (3.52)for the two-dimensional steady-state Advection-Diffusion equation were presented in Figure 4.1—Figure 4.11 and Table 4.7 respectively. The program was executed using the algorithm developed in MATLAB 7.8.0(R2009a) and a high speed computer hp Pavilion dv6 series of specifications: Processor—Intel (R) Core(TM) i5 CPU M520 @ 2.40 GHz, Memory (RAM)—4024 MB, 64 Bits Operating System, 240 GB Hard Disk, running Microsoft Window 7 Home Premium Operating system. After coding, the program was debugged, compiled and run. 3.6.2. COMPUTATIONAL ERRORS The main rationale of the numerical positive scheme of the Crank Nicholson approach was to compare the finite difference method to the results of the analytical solution of the Gaussian Plume Model. This was to show the level of accuracy and the deviation of this work from the exact solution. 𝐴𝐵𝑆𝐹𝐷𝑀𝑆𝐸𝑅𝑅𝑂𝑅 = 𝐴𝐵𝑆 𝐴𝑆𝑂𝐿 − 𝐹𝐷𝑀𝑆 𝐴𝑉𝐸𝑅𝐴𝐺𝐸𝐸𝑅𝑅𝑂𝑅𝑂𝐹𝐹𝐷𝑀𝑆 = 𝐴𝑉𝐸𝑅𝐴𝐺𝐸𝐹𝐷𝑀𝑆𝐸𝑅𝑅𝑂𝑅𝑆 𝐴𝑉𝐸𝑅𝐴𝐺𝐸𝐸𝑅𝑅𝑂𝑅𝑂𝐹𝐴𝑆𝑂𝐿 = 𝐴𝑉𝐸𝑅𝐴𝐺𝐸𝐴𝑆𝑂𝐿𝐸𝑅𝑅𝑂𝑅𝑆 𝐸𝑅𝑅𝑂𝑅𝐴𝑆𝑂𝐿/𝐹𝐷𝑀𝑆 = 𝐴𝑆𝑂𝐿 − 𝐹𝐷𝑀𝑆 These errors were computed using the Microsoft Excel and the output was displayed in Table 4.7. University of Ghana http://ugspace.ug.edu.gh 50 3.6.3. CALIBRATION, VALIDATION AND VERIFICATION The atmospheric dispersion model of Cs-137 was calibrated and validated by the Gaussian Plume analytical solution calculations and verified by comparing with results obtained from the Numerical solution. The quantity of mass concentration deposited at each Receptor zone driven by the wind was validated using input data from the safety analysis report and verified by comparison with the numerical results. University of Ghana http://ugspace.ug.edu.gh 51 CHAPTER FOUR RESULTS AND DISCUSSIONS The numerical and analytical solutions of Equation 3.8 (the advection- diffusion equation) subject to the boundary conditions and initial conditions were simulated and illustrated graphically. The influences of all the embedded parameters to the dispersion of radionuclide transport were considered. The parameters considered in the model were the drifting velocity (u), dynamic viscosity of air (𝜇), the emission rate (𝑄) and height (𝐻). 4.1. RESULTS The two accidental scenarios considered in this part were; when the source was at the ground level, (𝐻 = 0) and when the source was elevated above the ground, ( 𝐻 = 5𝑚). The plots in Figure 4.1 illustrated the maximum value of Cs-137 concentration at ground-level in both vertical and horizontal cross sections. The dispersion also showed that the Cesium-137 was swept from the maximum value of concentration into an elongated plume shape. The plot in Figure 4.1(a) depicted the dispersion of Cs-137 concentration in the plane when 𝑦 = 0, and demonstrated that the peak ground-level concentration occurs at the origin when 𝐻 = 0. The plot in Figure 4.1(b) illustrated the dispersion of Cs-137 in the horizontal cross-section(𝑧 = 0). The location of the peak ground-level concentration was marked with a red circle. University of Ghana http://ugspace.ug.edu.gh 52 4.1.1. SOURCE EMISSION AT GROUND–LEVEL (a) Vertical Cross-Section (𝑦 = 0) (b)Horizontal Cross-Section (𝑧 = 0) Figure 4.1: Cs-137 Dispersion Diagram at Ground-Level: (a)— Vertical Cross- Section (𝒚 = 𝟎,𝒖 = 𝟓 𝒎/𝒔); (b)— Horizontal Cross-Section (𝒛 = 𝟎, 𝒖 = 𝟓𝒎/𝒔). University of Ghana http://ugspace.ug.edu.gh 53 4.1.2. CONCENTRATION OF CESIUM-137 AT ELEVATED SOURCE Figures 4.2—4.3 illustrated the plots of Cs-137 concentration for an elevated source of height 5𝑚 at 5𝑚/𝑠. It was observed from the diagrams that when the source was elevated ( 𝐻 > 0), the concentration attains a maximum value of 𝐶𝑚𝑎𝑥 = 2𝑄/(𝜋𝑢𝐻 2) in Equation (3.33) along the plume centerline when 𝑦 = 0. The plots also showed that the mass concentration of Cesium-137 from a source height of 5𝑚 and wind speed, 5 𝑚/𝑠, was 4.3154 × 10−4 mg/m3. The graphs in Figure 4.2— 4.3 showed that Cs-137 concentration from the source decreases as the source height was elevated but increases at the Receptor zones as the source height was lowered. The steady state Gaussian equation (which was used in the implementation of the concentration) showed that the concentration was inversely proportional to the wind speed (the lower the wind speed the higher the concentration). University of Ghana http://ugspace.ug.edu.gh 54 Vertical Cross-Section (𝑦 = 0,𝑢 = 5𝑚/𝑠 ,𝐻 = 5m) Figure 4.2: Cs-137 Dispersion Diagram for Elevated Height: Vertical Cross-section (𝒛 = 𝟓𝐦, 𝒚 = 𝟎, 𝒖 = 𝟓 𝒎/𝒔). Map Legend Horizontal Cross-Section (𝑧 = 0, 𝑢 = 5𝑚/𝑠 ,𝐻𝑒𝑖𝑔ht = 5m) University of Ghana http://ugspace.ug.edu.gh 55 Figure 4.3: Cs-137 Dispersion Diagram for Elevated Height: Horizontal Cross-section (𝒚 = 𝟓𝒎, 𝒛 = 𝟎,𝒖 = 𝟓 𝒎/𝒔). Map Legend University of Ghana http://ugspace.ug.edu.gh 56 4.1.3. CESIUM-137 DISPERSION AT A SET OF RECEPTOR LOCATIONS 4.1.3.1. Numerical Solutions with Non-Zero Drift Velocity (𝒖 = 𝟐.𝟓𝒎/𝒔) Figure 4.4: Dispersion of Cs-137 Concentration at Multiple Receptor Zones: (𝒖 = 𝟐.𝟓𝒎/𝒔, 𝑯 = 𝟓𝐦) Figure 4.5: Mass of Cs-137 (mg) Deposited at each Receptor Zone: (𝒖 = 𝟐.𝟓𝒎/𝒔, 𝑯 = 𝟓𝐦) University of Ghana http://ugspace.ug.edu.gh 57 Figures 4.4 and 4.5 illustrated the dispersion of Cs-137 when the wind speed was 2.5m/s and Table 4.2 represented the mass of Cs-137 discharged at each receptor zone. Onyansana and Pokuase did not receive any Cs-137 concentration and the maximum concentration was deposited around Akokome Township. Table 4.1: Mass of Cs-137 Deposited at each Receptor Zone RECEPTORS Onyansana Pokuase Mpehuasem Nkotunsi Akokome Adenta Abokobi Teiman Mass (mg) 0 0 0.0004 0.0002 0.4686 0.0856 0.0580 0.0065 Total mass(mg) 𝟔.𝟏𝟗𝟎𝟏× 𝟏𝟎−𝟓 Map Legend North East The predominant wind direction in Accra is from South-West to North-East sectors. University of Ghana http://ugspace.ug.edu.gh 58 4.1.3.2. Numerical Solutions with Non-Zero Drift Velocity (𝒖 = 𝟓𝒎/𝒔) Figure 4.6: Dispersion of Cs-137 Concentration at Multiple Receptor Zones: (𝒖 = 𝟓𝒎/𝒔 ,𝑯 = 𝟓 𝒎) Figure 4.7: Mass of Cs-137 (mg) Deposited at each Receptor Zone: (𝒖 = 𝟕.𝟓𝒎/𝒔, 𝑯 = 𝟓𝐦) University of Ghana http://ugspace.ug.edu.gh 59 Table 4.2: Mass of Cs-137 Deposited at each Receptor Zone RECEPTORS Onyansana Pokuase Mpehuasem Nkotunsi Akokome Adenta Abokobi Teiman Mass (mg) 0 0 0.0012 0.0011 0.1487 0.0447 0.01103 0.0131 Total mass(mg) 𝟒.𝟐𝟗𝟒𝟏× 𝟏𝟎−𝟓 Figures 4.6 and 4.7 showed the dispersion of Cs-137 when the wind speed was 5m/s and Table 4.3 illustrated the mass of Cs-137 discharged at each receptor zone. The maximum mass concentration of 0.1487 mg was recorded at Akokome. University of Ghana http://ugspace.ug.edu.gh 60 4.1.3.3. Numerical Solutions with Non-Zero Drift Velocity (𝒖 = 𝟕.𝟓𝒎/𝒔) Figure 4.8: Dispersion of Cs-137 Concentration at Multiple Receptor Zones: (𝒖 = 𝟕.𝟓𝒎/𝒔, 𝑯 = 𝟓𝐦) Figure 4.9: Mass of Cs-137 (mg) Deposited at each Receptor Zone: (𝒖 = 𝟕.𝟓𝒎/𝒔, 𝑯 = 𝟓𝐦) University of Ghana http://ugspace.ug.edu.gh 61 Table 4.3: Mass of Cs-137 Deposited at each Receptor Zone RECEPTORS Onyansana Pokuase Mpehuasem Nkotunsi Akokome Adenta Abokobi Teiman Mass× 𝟏𝟎−𝟓 (mg) 0 0 0.0002 0.0001 0.9800 0.1500 0.1300 0.1340 Total mass(mg) 𝟏.𝟑𝟗𝟒𝟑× 𝟏𝟎−𝟓 Figures 4.8 and 4.9 showed the dispersion of Cs-137 when the wind speed was 7.5m/s and Table 4.4 illustrated the mass of Cs-137 discharged at each receptor zone. The least concentration, 0.0001 was recorded at Nkotunsi. The total concentration deposited at Onyansana and Pokuase was zero. This was due to the drift velocity of the wind from Southwest to the North-Eastern region. University of Ghana http://ugspace.ug.edu.gh 62 4.1.4. DISPERSION OF CESIUM-137 AT A SET OF RECEPTOR ZONES UNDER EXTREME CONDITIONS 4.4.1. Numerical Solutions under Extreme Condition (𝒖 = 𝟏𝟎𝒎/𝒔,𝑯 = 𝟓𝒎) Figure 4.10: Dispersion of Cs-137 Concentration: (𝒖 = 𝟏𝟎𝒎/𝒔, 𝑯 = 𝟓𝒎)) Figure 4.11: Mass of Cs-137 (mg) Deposited at each Receptor Zone: 𝒖 = 𝟏𝟎𝒎/𝒔 University of Ghana http://ugspace.ug.edu.gh 63 Table 4.4: Mass of Cs-137 Deposited at each Receptor Zone RECEPTORS Onyansana Pokuase Mpehuasem Nkotunsi Akokome Adenta Abokobi Teiman Mass(mg) 0 0 0.1361 0.0000 0.2935 0.6706 0.4815 0.0110 Total mass(mg) 𝟏.𝟓𝟗𝟐𝟖 × 𝟏𝟎−𝟓 Figures 4.10 –Figure 4.11 represented accident scenarios during extreme conditions where the wind speed was 10𝑚/𝑠.Under severe conditions, the Cs-137 concentration affected more receptor levels than during average situations. Moreover, it was observed that the Cs-137 concentration spread wider during these extreme conditions. 4.1.5. Analytical and Numerical Integration of Cs-137 Concentrations There is a correlation between the stability and the wind speed for the dispersion of contaminants. Similarly, the source height and the distance of the receptors can influence the discharge of the Cs-137 at the zones [38]. Table 4.6 and Figure 4.12 showed the correlation between these stabilities and variations. Table 4.6 showed the results obtained for the numerical and exact solutions. These established the association between concentration and stability indicated at Table 2.1, that stability played a major role in the dispersion of the Cs-137 concentration. Table 2.1 showed that, stability was slightly unstable at 5 m/s, a distance of 10 m, and this stability was confirmed in Table 4.6. The Graphs in Figure 4.12 showed that there was strong correlation between the analytical and the numerical solutions at the Neutral point in the Northeast region. University of Ghana http://ugspace.ug.edu.gh 64 Table 4.5: The comparison between the numerical and analytical concentrations Run no. Stability u(m/s) h(m) Distance × 𝟏𝟎−𝟐 (m) C( [mg/m3] Numerical Analytical 1 Very unstable (A) 1.0 0 500.0 4.6992× 10−5 4.8183 × 10−6 1 Very unstable (A) 1.5 0 500.0 3.1988 × 10−6 3.2122× 10−6 2 Slightly unstable (C) 2.0 1.0 500.0 2.3481 × 10−6 2.4055 × 10−6 2 Slightly unstable (C) 2.5 1.0 500.0 1.8385× 10−6 1.9244× 10−6 3 Moderately unstable (B) 3.0 2.0 500.0 1.5830× 10−6 1.5965 × 10−6 3 Moderately unstable(B) 3.5 2.0 500.0 1.3425× 10−6 1.3685 × 10−6 4 Slightly unstable C) 4.0 3.0 500.0 1.1743× 10−6 1.1885 × 10−6 4 Slightly unstable (C) 4.5 3.0 500.0 1.0372× 10−6 1.0564 × 10−6 5 Slightly unstable (C) 5.0 4.0 500.0 9.3640× 10−7 9.4091 × 10−7 5 Slightly unstable C) 5.5 4.0 500.0 8.5036× 10−7 8.5537× 10−7 6 Slightly unstable (C) 6.0 5.0 500.0 7.6518× 10−6 7.7362× 10−7 6 Slightly unstable (C) 6.5 5.0 500.0 7.1246× 10−7 7.1411 × 10−7 7 Neutral (D) 7.0 6.0 500.0 6.5107× 10−7 6.5231 × 10−7 7 Neutral (D) 7.5 6.0 500.0 6.0941× 10−7 6.0882 × 10−7 8 Neutral(D) 8.0 7.0 500.0 5.0044× 10−7 5.0044× 10−7 8 Slightly unstable C) 8.5 7.0 500.0 5.2033× 10−7 5.2687 × 10−7 9 Slightly unstable (C) 9.0 8.0 500.0 4.7462× 10−7 4.8658× 10−7 9 Moderately unstable (B) 9.5 8.0 500.0 4.6017× 10−7 4.6097 × 10−7 10 Moderately unstable (B) 10.0 9.0 500.0 4.1566× 10−7 4.2694 × 10−7 10 Moderately unstable (B) 10.5 9.0 500.0 4.0539× 10−7 4.0661 × 10−7 University of Ghana http://ugspace.ug.edu.gh 65 Figure 4.12: A Graph of Numerical Versus Analytical Solutions University of Ghana http://ugspace.ug.edu.gh 66 4.1.6. ERROR ANALYSIS The average errors for the two dimensional analysis of the Advection-Diffusion Equation for both the Analytical and the Numerical solutions were displayed in Table 4.7. Table 4.6: Error Analysis for ASOL & NSOL The errors encountered in this thesis were truncation errors, computational errors and round-off errors of parameter values. These errors caused the deviation of the analytical solutions from the numerical simulations as shown in Table 4.7. The average error was found to be 2.778094 %. Run No. ASOL NSOL ERROR %ERROR 1 4.82E-05 4.70E-05 1.19E-06 2.478052 2 3.21E-06 3.20E-06 1.21E-08 0.376689 3 2.41E-06 2.35E-06 5.74E-08 2.386198 4 1.92E-06 1.84E-06 8.59E-08 4.463729 5 1.60E-06 1.58E-06 1.91E-08 1.196367 6 1.37E-06 1.34E-06 2.60E-08 1.89989 7 1.19E-06 1.17E-06 2.22E-08 1.867901 8 1.06E-06 1.04E-06 2.02E-08 1.912154 9 9.41E-07 9.36E-07 4.51E-09 0.479323 10 8.55E-07 8.50E-07 5.01E-09 0.585711 11 7.74E-07 7.65E-07 9.11E-09 1.177581 12 7.14E-07 7.12E-07 2.01E-09 0.281469 13 6.52E-07 6.51E-07 1.21E-09 0.185495 14 6.09E-07 6.09E-07 0.00E+00 0.140122 15 5.00E-07 5.00E-07 4.40E-10 0.087923 16 5.27E-07 5.20E-07 6.56E-09 1.245089 17 4.87E-07 4.75E-07 1.20E-08 2.457972 18 4.61E-07 4.60E-07 8.00E-10 0.173547 19 4.27E-07 4.16E-07 1.13E-08 2.642057 20 4.07E-07 4.05E-07 1.22E-09 0.300042 AVE. ERR 4.86E-05 4.74E-05 1.20E-06 2.778094 University of Ghana http://ugspace.ug.edu.gh 67 4.2. DISCUSSION The source height influenced the dispersion of the Cesium concentration. The graphs (Figures 4.1—4.3) showed that as the source height increased, Cs-137 concentration diminished around the receptor zones. The largest ground-level concentration of Cs-137 was deposited at Akokome in the Northeast region from the source, with weaker peaks occurring at Nkotunsi. The Akokome Township was in the predominant wind direction (North-East wards) (Figure 4.4) from the source; hence the largest Cs-137 concentration was recorded in that vicinity. There was no Cs-137 concentration recorded at Onyansana and Pokuase and this was due to the location of these receptors and the prevailing wind from Southwest to North-Eastern region (Figure 4.4—4.6 ). The rate at which the radionuclide traveled to the receptors was dependent on the distance of the receptors from the released source and the wind velocity. The graphs (Figure 4.2—4.9) showed that the receptors which were close to the emitting source and were in the North-Eastern region of the wind received more concentration of Cs-137. Table 4.2 - Table 4.4depicted the Cesium-137 concentration at each receptor zone when the wind speed, 𝑢 was between 2.5 𝑚/𝑠 𝑎𝑛𝑑 7 𝑚/𝑠. From the diagrams, the values showed that, Onyansana, Pokuase and Nkotunsi, recorded no Cs-137 concentration from the source. This was due to the placement of the emitting source and the location of Onyansana, Pokuase and Nkotunsi. As the wind remained in one direction and propelled in the direction of the prevailing winds (SW to NE), then the Gaussian plume solution returns a value that is identically zero for any value of 𝑥 ≤ 0 (that is for any point lying in the North-West and South-West of the source. University of Ghana http://ugspace.ug.edu.gh 68 Figures 4.4—4.8 showed that Onyansana, Pokuase and Nkotunsi laid in the North- Western and South-Western sector from the source where the Gaussian plume solution was essentially zero at that point (Figure 2.2). The average activity concentration of Cs-137 was very low at Mpehuasem as compared to other receptor zones. This was due to the location of the receptor and the prevailing North- Easterly wind. The graphs in Figure 4.4 to Figure 4.9 showed that the maximum distance that was covered by Cesium-137 from the source was 1795m, 1621 m and 1550 at 2.5m/s, 5 m/s and 7 m/s respectively. Consequently, regulators must decontaminate as far as 4000m square meters Northeast wards during nuclear accident. Furthermore, since Akokome, Adenta, and Abokobi and in some cases Teiman recorded the maximum contaminant concentration levels respectively, the supervisory body should give much consideration to these receptor zones. Also, the receptors (Akokokme, Adentan, Abokobi, Teiman) that were located in the North-Eastern sector of the source normally received the highest contaminant concentration of Cs-137 as the source height reduced and at low wind speed. However, the weaker contaminant concentration was discharged around the receptors (Onyansana, Pokuase, Nkotunsi) which were located at the North-Western and South- Western part of the source respectively. Therefore, regulators must ensure that those receptors such as: factories, suburban communities and water bodies from the source should be protected and if possible relocate these receptors. In general, the ground- level concentration increased as the distance within the source increases, and then decreased rapidly beyond the source and as the wind speed increases, the ground level concentration decreases. University of Ghana http://ugspace.ug.edu.gh 69 CHAPTER FIVE CONCLUSIONS AND RECOMMENDATIONS 5.1. CONCLUSIONS The Gaussian Plume Model based on two dimensions Advection-Diffusion Equation had been developed. The model had three external parameters, namely Cs-137 diffusion coefficient, the drift velocity of air and the eddy diffusivity coefficient. The polluting source was characterized by location, emission rate and the period during which the source was active. The numerical method used to solve the Advection-Diffusion Equation model was the finite difference Crank-Nicholson method, and the algorithm was implemented using MATLAB 7.8.0(R2009a). From the study, it was observed that; (i) If the wind remained in one direction and blows in the direction of the prevailing North-East winds, then the Gaussian plume solution returns a value that is identically zero for any value of 𝑥 ≤ 0 ( that is for any point lying in the Southwest wards of the source); (ii) The ground-level concentration increased as the atmosphere becomes more stable (0.0012 mg/m3 − 0.1487 mg/m3); (iii) Reducing source height (ranging 5m — 2m) led to increasing in ground level concentration of Cs-137; (iv) Increasing wind speed (1m/s—5 m/s) caused Cs-137 concentration to decrease near ground level; (v) The drift velocity (in the range 5 × 10−4 to1 × 10−4m/s) of air moves Cs-137 from one region to another at the rate which is proportional to its magnitude. University of Ghana http://ugspace.ug.edu.gh 70 5.2. RECOMMENDATIONS (1). The volatile or semi-volatile gases expected to leak from the core during an accident can also be modeled in the same manner used for the present study. This could provide additional insight to the total contribution to external radiation due to each isotope, and can be used as a type of a screening and sensitivity analysis. The isotopes mentioned could include some of the noble radioisotopes of xenon and krypton and other halogens. (2). The study assumed that the vertical movement of Cs-137 has no significant effect, however in real situation for some cases this should be put in consideration (i.e. three dimensional models). University of Ghana http://ugspace.ug.edu.gh 71 REFERENCES [1]. Venkatram A., Introduction to AERMOD Fundamentals: Micrometeorology and Dispersion notes, Department of Mechanical Engineering, 2011, Retrieved March 12, 2013 from University of California on the World Wide Web, http:/www.engr.ucr.edu/~venky.com.. [2]. Cooper D. and Alley C., Air Pollution Control: The Design Approach,3rd edition, Waveland Press, Inc., Arizona., United State of America, 2011, pp. 56—105. [3] El-Harbawi M., Mustapha S. and Rashid A., Air Pollution Modeling, Simulationand Computational Methods, In Proc. A Review Paper presented at the ICERT: International Conference on Environmental Research and Technology, Penang, Malaysia, 3-5 May, 2007. [4]. Vermal T. and Ghosh A.K., Asian Nuclear Prospect 2010, Safety Division, Safety & Environment Group, Bhabha Atomic Research Centre, Mumbai-400, India, 2010. [5]. Akaho E., Maakuu B., Anim S., Emi-Reynolds G., BoaduH., Osae E., Akoto B., and Dodoo D., Ghana Research Reactor-1 Safety Analysis Report, GAEC, Kwabenya, Ghana, February 2006. [6]. Yang Yue-Wen, MNSR Training Manual – GHARR-1, Guangzhou, China, 1989. [7]. Fort G., Thyou H., Air Pollution and Particle Emissions: Nuclear Sciences and Applications, Research Journal of American Nuclear Society. Vol. 88, Issue 7, 1999, pp. 41-88. University of Ghana http://ugspace.ug.edu.gh 72 [8]. Environmental Protection Agency., Air pollution, Retrieved on Wednesday, 6th July, 2012. http://www.epa.gov/ebtpages/airairpollutants.html.com. [9]. Tuner D.B., Workbook of Atmospheric Dispersion Estimates: An introduction to Dispersion Modeling, 2nd edition, CRC press, USA, 2001, pp. 32-56. [10]. http://en.wikipedia.org/wiki/three_mile_island. Retrieved on Monday, 16th May, 2012. [11]. http://en.wikipedia.org/wiki/Fukushima_Daiichi_nuclear_disaster. Retrieved on Monday, 16th May, 2012. [12]. Chakraborty M., Ahmad M., Singh R., Pal D., and Bandopadhyay C. K., Determination of the emission rate from various opencast mining operations; Environmental Modeling and Software, Vol.3, No.7.,2009, pp. 467—480, 705. [13]. McRae D., Gregory J., and William R., Numerical Solution of the Atmospheric Diffusion Equation for Chemically Reacting Flows, Journal of Computational Physics, 2001, pp.338-454. [14]. Arkhipov B., Vladimir K., Viacheslav S., Dmitry S. and Yulia Y., Numerical Modeling of Pollutant Dispersion and Oil Spreading by the Stochastic Discrete Particles Method Studies in Applied Mathematics, 2010, pp. 87-104. [15]. Noel H. and Stephen P., Stabilization of the Evolutionary Convection Diffusion Problem: Introduction and Experiments, Millersville University of Pennsylvania, USA, 2007, pp. 115—254. [16]. Osalu A., Fatehifar E., Kaynejad M. and Elkamel A., Developing a new air pollution dispersion model with chemical reactions based on multiple cell approach, Second International Conference on Env. and Computer Science, Paper presented at the 4th ULCOS seminar, , Germany, 1-2 October, 2008. University of Ghana http://ugspace.ug.edu.gh 73 [17]. Reeden R., Significant Dust Dispersion Models for Mining Operations, Report for National Institute for Occupational Safety and Health, Research Laboratory, Pittsburgh, USA., IC 9478 information circular, 2005. [18]. Venkatram S. and Seinfeld A., A Transport Pollution Investigation in Vilnius Streets, Journal of Environmental Engineering and Landscape Management, Vol.42, No. 13, 2011, pp.54-59. [19]. McKibbin R., Lim L., Smith T. and Winston L., A Model for Dispersal of Eruption Ejecta: Proceedings World Geothermal Congress Massey University, Albany Campus, Private Bag 102 904, NSMC, Auckland, New Zealand, 2005, pp.205-310. [20]. Schonar J., Environmental Modeling: Fate and Transport of Pollutants in Water Air and Soil, John Wiley & sons, Inc., New York, 2007. [21]. Indra N., Pascal. and Gurdeep S., Classification of Air Pollution Dispersion Models: Proceedings of the National Seminar on Environmental Engineering with special emphasis on Mining Environment, Journal of Environmental Engineering and Landscape Management, New York 11(2), 2004, pp.54-59. [22]. Perry S.G., A Dispersion Model for Industrial Source Applications. Part II: Model Performance against 17 Field Study Databases; Applied Meteorology, 2009, pp.695-707. [23]. Burden G., Richard L. and Douglas F., Numerical Analysis, Boston, Massachusetts: PWS-KENT publishing company, 2008. [24]. KanevceG.,Dispersion Modeling for Regulatory Applications Thermal Science, Vol. 10, No. 6, 2006, pp. 141-154. [25]. Draxler R., Across North America Tracer Experiment (ANTEX): Sampling and Analysis, Atmospheric Environment, Boston, Massachusetts, USA, 2006. University of Ghana http://ugspace.ug.edu.gh 74 [26]. Greenhill R., Particle Dispersion Models for Mining Operations, 2nd edition, Delton Press, Cambridge, Massachusetts, USA, 2003, pp. 45—78. [27]. Mihaiella C., Victoria T., Silviu I. and Adina I., Case Study on Pollution Prediction through Atmospheric Dispersion Modeling: US AEC Report ANL- 5673, Acoustics and Air Quality Laboratory, 2010. [28]. Colberga C., Tonab W., Stahelb M., Meierc J. and Staehelina B., Comparison of road traffic emission model with emissions derived from measurements in the Gubrist road tunnel, International Journal of Mechanical Science, Vol. 45. No 8, Switzerland, Atmospheric Environment, 2005, pp. 39—112. [29]. Bierly, E. W. and Hewson E. W., Some restrictive meteorological conditions to be considered in the design of nuclear plants, Journal of Applied Meteorology, 1990, pp.38-89. [30]. Nappo H., and Williams D., Twenty years experience with post-Chernobyl Thyroid Cancer, Edition 5, The Routney and Hill Company, Barberton, Ohio, 2007, pp.65-102 [31]. Wikipedia, Pollutant and greenhouse gases, Retrieved on Wednesday, 6thJuly, 2011, http://en.wikipedia.org/wiki/pollutant_green_house_gases. [32]. MichealT., Scientific Computing: An Introductory Survey, 2ndedition, The McGraw-Hill Companies, Inc., ISBN, Washington D.C., 2005. [33]. Vannucci M., Coolla V. and Haines N., Air Dispersion Modeling for the Assessment of ULCOS Technologies, Paper presented at the 4th ULCOS seminar, France, 1-2 October 2008. [34]. Turner D., Bruce M.,Workbook of Atmospheric Dispersion Estimates, 2nd edition, Delton Press, Cambridge, Massachusetts, USA, 2003, pp. 45—78. University of Ghana http://ugspace.ug.edu.gh 75 [35]. Recktenwald G., Finite-Difference Approximations to the Heat Equation, Portland State University, Portland, Oregon, 2004. [36]. Lanouette J.,Kemeny Commission Report, Bulletin of the Atomic Scientists, Vol. 36, Issue 1, 1992, p. 45. [37]. Chock A., David P., and Dunker M.K., A comparison of numerical methods for solving the advection equation, Atmospheric Environment, 1996, pp.11- 24. [38]. StockieJ., The Mathematics of Atmospheric Dispersion Modeling; Vol. 1, No. 2, 2007. [39]. William F., Numerical Methods for Partial Differential Equations, Academic Press, Inc., 3rd Edition, Boston, 1992. [40]. Winterton R.S., Thermal Design of Nuclear Reactors, Pergamon Press, New York, USA, 1981, pp. 110-250. [41]. Dop J., Han K., and Van G., Buoyant Plume Rise in a Lagrangian Framework, Atmospheric Environment, Vol. 7, No. 8, 2008, pp. 35-134. University of Ghana http://ugspace.ug.edu.gh 76 APPENDICES APPENDIX I TABLE OF LAPLACE TRANSFORMS 𝒇(𝒕) 𝒇 (𝒔) 1. 𝛿(𝑡) 1 2. sinh⁡(𝑎𝑡) 𝑎 𝑠2 − 𝑎2 3. cosh⁡(𝑎𝑡) 𝑎 𝑠2 − 𝑎2 4. 1 𝜋𝑡 𝑒−𝑎 2/4𝑡 1 𝑠 𝑒−𝑎/𝑠 (𝑎 ≥ 0) 5. 𝑓 ′(𝑡) 𝑠𝑓 (𝑠)− 𝑓 (0) 6. 𝑓 ′′ (𝑡) 𝑠2𝑓 𝑠 − 𝑠𝑓 0 − 𝑓 ′(07) 7. 𝑓 𝑡 − 𝑎 ℋ(𝑡 − 𝑎) 𝑒−𝑎𝑠𝑓 𝑠 (𝑎 ≥ 0) 8. 𝑒𝑎𝑡𝑓(𝑡) 𝑓 (𝑠 − 𝑎) 9. 𝑒−𝑎𝑡𝑓(𝑡) 1 𝑠 + 𝑎 University of Ghana http://ugspace.ug.edu.gh 77 𝑐 𝑥, 𝑦, 𝑧 = 𝑄 𝑢 𝑎 𝑥,𝑦 . 𝑏(𝑥, 𝑧) For 𝑎 𝑥,𝑦 the Laplace transform in 𝑟 is: 𝜌𝑎 − 𝑎 0,𝑦 = 𝜕2𝑎 𝜕𝑦2 Laplace transform in 𝑦 is: 𝜂2𝑎 − 𝜂𝑎 𝜌, 0 − 𝜕𝑎 𝜕𝑦 𝜌, 0 − 𝜌𝑎 = −1 Inverse transform in 𝜂 to obtain: 𝑎 𝜌, 𝑦 = 𝐶1 cosh 𝜌𝑦 − 𝐶2 2 𝜌 𝑒 𝜌𝑦 − 𝑒− 𝜌𝑦 Inverse transform in 𝜌 to get 𝑎(𝑥,𝑦) = 𝐶2 𝜋𝑥 𝑒𝑥𝑝 −𝑦2 2𝑥 𝑎 𝑟,𝑦 = 1 2𝜋𝑥 𝑒−𝑦 2/2𝑥 For 𝑏(𝑟, 𝑧) the Laplace transform of the separable form in 𝑟 is given by: 𝜕2𝑏 𝜕𝑧2 − 𝜌𝑏 = 𝛿(𝑧 − 𝐻) Laplace transform in: 𝑏 𝜌, 𝜁 = ℒ𝑧 𝑏 𝜌, 𝑧 University of Ghana http://ugspace.ug.edu.gh 78 Inverse transform in 𝜁 to obtain: 𝑏 𝜌, 𝑧 = 𝑏 𝜌, 0 𝑐𝑜𝑠𝑕 𝜌𝑧 − 1 𝜌 𝑠𝑖𝑛𝑕 𝜌(𝑧 − 𝐻) b ρ, z = 1 2 ρ e− ρ(z−H) + e− ρ(z+H) Inverse transform in 𝜌 yields: b x, z = 1 2πr e−(z−H) 2 2x + e−(z+H) 2 2x University of Ghana http://ugspace.ug.edu.gh 79 APPENDIX II %======================================================== % THIS MATLAB CODE, "DASIM" IS A SUBROUTINE TO THE MAIN PROGRAM THAT % SETS UP VARIOUS PHYSICAL PARAMETERS FOR THE ATMOSPHERIC DISPERSION % PROBLEM [38] %-------------------------------------------------------------------- %Source emission source data: source.n =1 ; % # of sources source.x = [250]; % x-location (m) source.y = [ 0]; % y-location (m) source.z = [ 0]; % height (m) source.label=['SO']; source.Q = [35] * tpy2kgps; % emission rate (kg/s) % Set locations of receptors where deposition measurements are made: recept.n = 8; % # of receptors recept.x = [ 60, 76, 66, 456, 514, 654, 2345, 110 ]; % x location (m) recept.y = [ 34, 78, 32, 128, 12, 22, 33, 25 ]; % y location (m) recept.z = [ 2, 5, 10, 31, 10, 11, 5, 11 ]; % height (m) recept.label=[ ' Nsanwam '; ' Pokuase '; ' Kwasima '; ' Achiasi '; ' Abokobi ';... 'Apenkwa '; ' Afianya '; ' Mampong ' ]; % Contaminant parameters (Ceasium-137): G = 9.8; % gravitational acceleration (m/s^2) D = 1.8e-5; % dynamic viscosity of air (kg/m.s) Den = 6000; % density of Cs-137 (kg/m^3) R = 0.80e-6; % diameter of Cs particles (m). Depo = 0.0092; % Cs-137 deposition velocity (m/s) Rset = 2*rho*grav*R^2 / (9*mu); % settling velocity (m/s): Mr= 65.4e-3; % molar mass of Cs (kg/mol) dia = 0.162; % receptor diameter (m) University of Ghana http://ugspace.ug.edu.gh 80 APPENDIX III %==================================================================== % The code "Conct-Mass" solves the atmospheric dispersion problem. . %-------------------------------------------------------------------- -- %(x=500m, y=150m at Source Height=5m, u=5m/s) DASIM; % read parameters from a file Uwind =10 ; % wind speed (m/s) % Set plotting parameters. dep = 0; warning( 'OFF', 'MATLAB:divideByZero' ); for i = 1 : source.n, % Sum up ground-level Cs concentrations from each source at all mesh % points, shifting the (x,y) coordinates so the source location is % at the origin. glc = glc + ermak( xmesh-source.x(i), ymesh-source.y(i), 0.0, ... source.z(i), source.Q(i), Uwind, Wdep, Wset ); px = 100; py = nx; xlim = [0, 2000]; ylim = [-100, 400]; x0 = xlim(1) + [0:nx]/nx * (xlim(2)-xlim(1)); % distance along wind direction (m) y0 = ylim(1) + [0:ny]/ny * (ylim(2)-ylim(1)); % cross-wind distance (m) [xmesh, ymesh] = meshgrid( x0, y0 ); % Do the same for the concentration at each receptor location and % scale the result by (A * dt * Wdep) to obtain a total deposition % (in mg) over the time interval 'dt'. dep = dep + (A * dt * Wdep) * ... ermak(recept.x-source.x(i), recept.y-source.y(i), recept.z, .. source.z(i), source.Q(i), Uwind, Wdep, Wset ); end warning( 'ON', 'MATLAB:divideByZero' ); % Plot contours of ground-level Cs concentration. figure(1) clist = [ 0.001, 0.01, 0.02, 0.05, 0.1 ]; glc2 = glc*1e6; % convert concentration to mg/m^3 University of Ghana http://ugspace.ug.edu.gh 81 [c2, h2] = contourf(xmesh, ymesh, glc2, clist ); % axis equal % (for plots in paper) clabel(c2, h2, 'FontSize', smallfont-2 ) colormap(1-winter) % These colors make the labels more readable %colorbar set(gca, 'XLim', xlim ), set(gca, 'YLim', ylim ) xlabel('x [miles(*10^ -2)]'), ylabel('Concentration(Kg/m^3)') title(['Cs-137 concentration (mg/m^3), max = ', sprintf('%5.2f', max(glc2(:)))]) grid on % Draw and label the source and receptor locations. hold on plot(source.x, source.y, 'ro', 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'r' ) text(source.x, source.y, source.label, 'FontSize', smallfont, 'FontWeight','bold' ); greencolor = [0,0.816,0]; plot(recept.x, recept.y, 'g^', 'MarkerEdgeColor', 'k', 'MarkerFaceColor', greencolor ) text(recept.x, recept.y, recept.label, 'FontSize', smallfont, 'FontWeight', 'bold' ) hold off shg print -djpeg 'glc.jpg' % Generate a second plot showing mass of Cs-137 deposited at each receptor. dep2 = dep * 1e6; % deposition in mg fprintf( 1, '\nDeposition in each town (mg):\n' ); disp(dep) fprintf( 1, 'Total deposited in all town zones (mg):\n' ); disp(sum(dep)) University of Ghana http://ugspace.ug.edu.gh 82 APPENDIX IV %======================================================== % THE CODE "G-LEVEL" PLOTS THE PLUME WITHOUT RECEPTOR ZONES. THE CONTAMINANT CONCENTRATION AND THE GROUND LEVELS OF THE Cs-137 %==================================================================== % Q = 10; % source emission rate (kg/s) U = 5; % wind speed (m/s) H = 1; % source height (m) K = 1; % eddy diffusion coefficient (m^2/s) rmax = 10; % size of domain in r (m) xmax = rmax*U/K; ymax = 0; % size of domain in y (m) zmax = 0; % size of domain in z (m) N = 400; % number of plotting points % Ground-level maximum concentration: x0 = U*H^2/(4*K) % x coordinate c0 = 2*Q/(pi*U*H^2*exp(1)) % concentration value % Plot the centerline concentration (y=0): [rr,zz] = meshgrid( [1:N]*rmax/N, [1:N]*zmax/N ); xx = rr*U/K; cc = Q./(4*pi*U*rr) .* ... ( exp(-(zz-H).^2./(4*rr)) + exp(-(zz+H).^2./(4*rr)) ); clist = [1e-5, 1e-4, 1e-3, 1e-2, 0.025, 0.05, 0.1]; [cs,ch] = contourf( xx, zz, cc, clist ); colormap jet clabel(cs) xlabel('x'), ylabel('z') hold on plot(x0+eps, 0+eps, 'ko') plot(0+eps, H+eps, 'rs') %arrow( [0.05*xmax, 0.8*zmax], [0.2*xmax,0.8*zmax], 'Width', 2 ) text( 0.05*xmax, 0.88*zmax, 'wind' ) hold off University of Ghana http://ugspace.ug.edu.gh 83 APPENDIX V %------------------------------------------------------------------- % This code, "Num-Ana.m" compares the crosswind integrated of the Cs-137 Concentration (mg/m^3). %------------------------------------------------------------------- nx=150; % geometry of the contaminant dispersion w=20; %(units on the x-axis) vel=5;% velocity of wind ( negative or positive) sigma=13.1; ampl=5.0; %(height of the source) nt=20;% Time step dt=1e-2; %Dispersion Rate dx=w/(nx-1); dx=w/(nx-1); x=0:dx:w; xc=21; % Plume Centreline t=ampl*exp(-(x-xc).^2/sigma^2); vx=ones(1,nx)*vel; foriteme=1:nt for ix=2:nx-1 %Tnew(ix)=0.1; end Tnew(1)=t(1); %Tnew(nx)=1; %t=Tnew(1); time=iteme*dt; T_anal=ampl*exp(-(x-xc-time*vel).^2/sigma^2); figure(1),clf,plot(x,t,x,T_anal), legend('analytical','numerical') xlabel('x(m)'), ylabel('Concentration(Kg/m^3)') drawnow end University of Ghana http://ugspace.ug.edu.gh 84 %==================================================================== % Calibration and Verification Matlab Codes %-------------------------------------------------------------------- % Input parameters: % %x=50;% - receptor locations: distance along the wind direction, with % the source at x=0 (m) %y=1;% - receptor locations: cross-wind direction (m) %z=1;% - receptor locations: vertical height (m) %H=5;% - source height (m) %Q=10;% - contaminant emission rate (mg/s) %U=8;% - wind velocity (m/s) %Wset=9.8;% - gravitational settling velocity (m/s) %Wdep=20;% - deposition velocity (m/s) % % First, define the cut-off velocity, below which concentration = 0. Umin = 0.0; sigmay = ay*abs(x).^by .* (x > 0); sigmaz = az*abs(x).^bz .* (x > 0); % Calculate the eddy diffusivity (m^2/s). Kz = 0.5*az*bz*U*abs(x).^(bz-1) .* (x > 0); % K = 0.5*U*d(sigma^2)/dx % Calculate the contaminant concentration (kg/m^3) using Ermak's formula. if U