University of Southern California Dissertations and Theses

Modeling of Earth's ionosphere and variational approach for data assimilation

(USC Thesis Other)

Modeling of Earth's ionosphere and variational approach for data assimilation

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content MODELING OF EARTH’S IONOSPHERE AND V ARIATIONAL APPROACH FOR DATA ASSIMILATION by Vardan Akopian A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (APPLIED MATHEMATICS) May 2008 Copyright 2008 Vardan Akopian Dedication To my parents ii Acknowledgements I would like to express my deep and sincere gratitude to my advisor, Dr. Chunming Wang. I could not have imagined having a better advisor and mentor for my disserta- tion. Always abundantly helpful, understanding, encouraging, Dr. Wang helped and guided me throughout this thesis, and gave me many lessons which will be invaluable in my future career. I would like to eagerly thank Brian Wilson who shared with me lots of his insights and helped in many stages of the GAIM++ development as well as my thesis. Our endless late evening discussions were huge stimuli for my mind and encouraged me to continue and finish this work. My sincere thanks are due to Tony Mannucci, supervisor of JPL’s Remote Sensing group, for giving me the opportunity to work at JPL, and thus have access to all the knowledge, data and the tools that were so helpful during this work. This work would not have been possible without the support of all my friends at the JPL’s Remote Sensing group, Tony Mannucci, Brian Wilson, Xiaoqing Pi, George Hajj, Attila Komjathy, Byron Iijima, Lukas Mandrake. Their wast knowledge in the domain of the ionospheric research was an endless source of learning for me. One last thank you goes to my family and friends who always supported me during these years, and especially to my parents who have ignited my love for science and implanted in me the believe that science matters in all times, good or bad. iii Table of Contents Dedication ii Acknowledgements iii List of Figures vi Abstract ix Chapter 1: Introduction 1 Chapter 2: Plasma physics of the ionosphere 6 2.1 Basic equations of an ionosphere model . . . . . . . . . . . . . . . . 7 2.1.1 The law of balance of mass . . . . . . . . . . . . . . . . . . . 7 2.1.2 The law of balance of momentum . . . . . . . . . . . . . . . 11 2.2 Complete equation for ionospheric dynamics in differential form . . . 17 2.3 Models of the Earth’s geomagnetic field . . . . . . . . . . . . . . . . 18 2.4 Partial differential equations for ionospheric dynamics . . . . . . . . 21 Chapter 3: Well posedness of equations 26 3.1 Existence and uniqueness of the solution - Single ion case . . . . . . . 29 3.2 Existence and uniqueness of the solution - multiple ion case . . . . . 34 Chapter 4: Formulation of numerical approximation 35 4.1 Space discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.2 Numerical form of the equations . . . . . . . . . . . . . . . . . . . . 41 4.3 Time discretization and Stability analysis . . . . . . . . . . . . . . . 49 Chapter 5: Validation of the forward model 53 5.1 Characteristics of the model ionosphere . . . . . . . . . . . . . . . . 54 5.2 The sensitivity of the model to some physics drivers . . . . . . . . . . 65 5.2.1 ExB Drift . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.2.2 Neutral Wind . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2.3 Solar Production . . . . . . . . . . . . . . . . . . . . . . . . 71 Chapter 6: The 4DV AR Process 73 6.1 The 4DV AR equations . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.2 4DV AR in GAIM++ . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6.2.1 Parametrizations for the drivers . . . . . . . . . . . . . . . . 78 6.2.2 4DV AR difficulties and issues . . . . . . . . . . . . . . . . . 83 iv Chapter 7: Validation of the 4DV AR Method 85 7.1 OSSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.1.1 Solar Production . . . . . . . . . . . . . . . . . . . . . . . . 88 7.1.2 ExB Drift . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 7.1.3 Neutral Wind . . . . . . . . . . . . . . . . . . . . . . . . . . 92 7.1.4 Solar Production and ExB Drift together . . . . . . . . . . . . 96 7.1.5 ExB Drift and Neutral Wind together . . . . . . . . . . . . . 97 Chapter 8: Design and Implementation of GAIM++ 104 8.1 Basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 8.1.1 Memory and other resource management . . . . . . . . . . . 104 8.1.2 Input/Output . . . . . . . . . . . . . . . . . . . . . . . . . . 105 8.1.3 Loader/Generator apparatus . . . . . . . . . . . . . . . . . . 106 8.1.4 Other utilities . . . . . . . . . . . . . . . . . . . . . . . . . . 109 8.2 Geometry Module . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 8.3 Physics Input Module . . . . . . . . . . . . . . . . . . . . . . . . . . 115 8.4 Forward Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 8.5 Observation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 8.6 4DV AR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 8.7 Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 Chapter 9: Conclusion 132 References 134 v List of Figures 1.1 Map of occultations locations in one day . . . . . . . . . . . . . . . . 3 1.2 GPS network map . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 p andq lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2 pq region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.1 Uniform grid inp andq . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2 Minimum, maximum and midlat values ofp . . . . . . . . . . . . . . 39 4.3 A low resolution GAIM++ grid . . . . . . . . . . . . . . . . . . . . . 40 4.4 GAIM++ grid lower boundary . . . . . . . . . . . . . . . . . . . . . 42 5.1 Ion density profiles, solar maximum year . . . . . . . . . . . . . . . . 55 5.2 Ion density profiles, solar minimum year . . . . . . . . . . . . . . . . 56 5.3 Ion density profiles from other sources . . . . . . . . . . . . . . . . . 57 5.4 O + density longitude slices . . . . . . . . . . . . . . . . . . . . . . . 59 5.5 H + density longitude slices . . . . . . . . . . . . . . . . . . . . . . . 60 5.6 He + density longitude slices . . . . . . . . . . . . . . . . . . . . . . 61 5.7 VTEC map, solar maximum year . . . . . . . . . . . . . . . . . . . . 62 5.8 VTEC map, solar minimum year . . . . . . . . . . . . . . . . . . . . 63 5.9 GIM VTEC map, solar maximum year . . . . . . . . . . . . . . . . . 63 5.10 GIM VTEC map, solar minimum year . . . . . . . . . . . . . . . . . 64 5.11 EB drift pre-reversal enhancement perturbation . . . . . . . . . . 66 5.12 Sensitivity ofO + toEB drift . . . . . . . . . . . . . . . . . . . . 66 5.13 Sensitivity of VTEC toEB drift . . . . . . . . . . . . . . . . . . 67 5.14 Wind perturbation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.15 Sensitivity ofO + to neutral wind . . . . . . . . . . . . . . . . . . . . 69 vi 5.16 Sensitivity of VTEC to neutral wind . . . . . . . . . . . . . . . . . . 70 5.17 Sensitivity ofO + to solar production . . . . . . . . . . . . . . . . . . 71 5.18 Sensitivity of VTEC to solar production . . . . . . . . . . . . . . . . 72 6.1 Cubic spline basis functions forEB drift perturbation . . . . . . . 82 7.1 Cost function, solar production . . . . . . . . . . . . . . . . . . . . . 89 7.2 Prefit and postfit residuals, solar production . . . . . . . . . . . . . . 89 7.3 Convergence, solar production . . . . . . . . . . . . . . . . . . . . . 90 7.4 Cost function,EB drift . . . . . . . . . . . . . . . . . . . . . . . 91 7.5 Prefit and postfit residuals,EB drift . . . . . . . . . . . . . . . . . 91 7.6 EB pattern convergence . . . . . . . . . . . . . . . . . . . . . . . 92 7.7 Cost function, neutral wind . . . . . . . . . . . . . . . . . . . . . . . 93 7.8 Prefit and postfit residuals, neutral wind . . . . . . . . . . . . . . . . 93 7.9 Neutral wind pattern convergence . . . . . . . . . . . . . . . . . . . 95 7.10 Cost function, solar production andEB drift together . . . . . . . 96 7.11 Prefit and postfit residuals, solar production andEB drift together . 96 7.12 Production andEB drift pattern convergence . . . . . . . . . . . . 97 7.13 Cost function,EB drift and neutral wind together . . . . . . . . . 98 7.14 Prefit and postfit residuals,EB drift and neutral wind together . . . 98 7.15 EB pattern at the end of cycle 8 . . . . . . . . . . . . . . . . . . . 99 7.16 Neutral wind pattern at the end of cycle 8 . . . . . . . . . . . . . . . 99 7.17 Vertical TEC maps for theEB and neutral wind estimation . . . . 100 8.1 MatrixFile class hierarchy . . . . . . . . . . . . . . . . . . . . . 105 8.2 The flow of the forward model . . . . . . . . . . . . . . . . . . . . . 117 8.3 Observation class hierarchy . . . . . . . . . . . . . . . . . . . . 119 8.4 The flow of one cycle of the 4DV AR process . . . . . . . . . . . . . . 120 vii 8.5 ThePerturbation class hierarchy . . . . . . . . . . . . . . . . . 123 8.6 The flow of the Kalman filter process . . . . . . . . . . . . . . . . . . 129 viii Abstract The Earth’s ionosphere plays an important role in wireless communication. One of the primary goals of this research is to improve the techniques for monitoring and forecasting of the conditions of the Earth’s ionosphere. The work in this thesis focuses on modeling of the ion and electron density changes in the Earth’s ionosphere and the assimilation of remote sensing measurements into the numerical model. The governing equations of the ionosphere consist of equations of fluid for chemically active plasma in the ionosphere. These equations form a system of coupled nonlinear convection- diffusion equations. In this thesis, we establish the well-posedness of the equations in some special cases. We also establish the numerical stability and consistency of the approximation. The resulting numerical model of the ionosphere produces stable and physically meaningful predictions of densities for all major ion species in the iono- sphere. The numerical method is used in this research as a basis for the development of a 4-dimensional variational (4DV AR) data assimilation approach to help in the esti- mation of the driving forces in the ionosphere. The 4DV AR technique uses the adjoint approach in the computation of the gradient of performance functional. In addition to theoretical discussions on the mathematical properties of the basic ionospheric model and its numerical implementation, extensive work has been devoted to the validation of both the basic ionospheric model and 4DV AR implementation. In this work we demonstrate that the new model can accurately reproduce all major char- acteristics of the ionosphere and the 4DV AR approach can be effective in determining the crucial drivers of the Earth’s ionosphere. ix Chapter 1 Introduction The Earth’s ionosphere is a layer of atmosphere at an altitude of 60km to 2000km above the Earth’s surface. The existence of ionosphere was discovered with the devel- opment of wireless communication through radio broadcast. The presence of ionized gases reflects radio waves to allow transmission over long distance in some wave- length. The impact of the Earth’s ionospheric conditions on social, economical and national defense activities has become increasingly significant as satellite navigation and wireless communication systems become part of our daily lives. For example, significant degradation of positioning accuracy in GPS-based navigation systems can be caused by disturbed ionospheric conditions. Similarly, the propagation of a High Frequency (HF) radio wave can be strongly affected by ionospheric irregularities. The need for a reliable and accurate ionospheric specification and forecast system is widely recognized by the space weather community. Indeed, the National Space Weather Pro- gram (NSWP) implementation plan drafted in 1995 has set the understanding of the day-to-day variability of the large-scale ionospheric features and small-scale plasma density irregularities as one of major goals of the NSWP. In 1999, the Department of Defense identified global ionospheric data assimilation as one of 12 topics for the Multidisciplinary University Research Initiative (MURI) program. Data assimilation is a process of reconciling a first principle based model of a complex system with actual observations of the system to obtain accurate estimate of the initial state and, in some case, relevant parameters of the system [13]. The early 1 development of data assimilation approach is primarily in the area of weather fore- cast. In this case, the system includes the atmospheric temperature, moisture, pressure and velocity as states. Since the mid-1990s, a collaborative research team of scientists from the Jet Propulsion Laboratory (JPL) and the University of Southern California has been using data assimilation techniques for ionospheric monitoring and forecast- ing. As one of two teams supported by the 1999 MURI effort, the USC/JPL team has developed an early version of a Global Assimilative Ionospheric Model (GAIM) [32], [14], [25], [21]. This model has demonstrated the ability to incorporate a diverse source of ionospheric measurements to significantly improve the accuracy of estima- tion of ionospheric state variables. In fact one of goals of the early development of the GAIM is to demonstrate that consistent information on ionosphere can be derived using vastly different measure- ments. These measurements include Total Electron Content (TEC) along line of sight between GPS satellites and ground based receivers; TEC along line of sight between a Low Earth Orbit (LEO) satellite and GPS satellites, also known as GPS occultation data; Vertical Total Electron Content (VTEC) from space based altimeters such as TOPEX and JASON; In-situ electron density measured by DMSP satellites; Vertical electron density profiles from ground based radar also known as ionosonde; 2 Earth limb profiles of ultra-violet airglow emissions measured from DMSP satel- lite instruments (SSULI), and Horizontal scan of ultra-violet airglow emissions measured from DMSP satellite instruments (SUSSI). As we can see, the majority of these measurements relies on space based instrument. Indeed the most important stimulant for ionospheric data assimilation is the rapid increase in the available measurements of the ionosphere. For an example, with launch of a family of 6 LEO satellites COSMIC in 2006, the number of GPS occultation mea- surements provides complete global cover of the ionosphere as shown in figure 1.1. For each point in this figure, a vertical profile of electron density can be retrieved using the empirical method of Abel inversion. A data assimilation model has the potential to make use of the same data without approximating assumptions. Another example of high fidelity and global coverage ionospheric measurement is the ground based TEC measurements. The JPL group has cooperative data exchange agreement with a net- work of GPS receiver station operators to collect global slant TEC measurements every 5 minutes. In figure 1.2, the coverage of the GPS receiver network is shown. Figure 1.1: One day coverage of 6000 occultations from Oersted, SUNSAT, SAC-C, CHAMP, GRACE, and COSMIC. Each occultation is represented by one dot 3 Figure 1.2: A map showing the distribution of selected stations of the global GPS network. The circles around the stations enclose projection areas at 400 km altitude, in which GPS data are taken at elevation angles higher than 10 degrees. The early version of the ionospheric model was intentionally simplified. In fact, it models only the most predominant ion speciesO+ in F-region of the ionosphere. As a prototype for GAIM, it was sufficient. However, important ionospheric physics is lack- ing in the early version of the model. The work of this thesis extends the early GAIM in developing a complete multi-ion version of the basic ionospheric model. The pri- mary challenge in this extension is that the multi-ion version of the ionosphere model equations is a nonlinear degenerated system of parabolic partial differential equation while the single ion version of the model equations is a scale linear partial differential equation. The method of approximating the solution of the nonlinear system of PDE requires the maintenance of the stability and positiveness of the solution. A second aspect of this thesis consists of further development of variational approach for data assimilation. Indeed the earlier experience with ionospheric data assimilation has used nearly exclusively recursive estimation techniques such as 4 Kalman filter. However, one most challenging aspect of ionospheric data assimila- tion is the estimation of the ionospheric driving forces such as the EB ion drift velocity, neutral wind velocity, neutral particle density and solar radiation intensity. We have experimented with the earlier version of variational data assimilation tech- niques (4DV AR). In this thesis a much more systematic exploration of the 4DV AR technique is carried out. An outline of the thesis is as follows. In chapter 2 we present the physiccal laws governing the ionosphere and the mathematical equations that model these laws. In chapter 3 we show the well-posedness of these equations under some special assump- tions. In chapter 4 we described the numerical approximation for the ionosphere equa- tions used in our implementation. Chapter 5 focuses on validating the results obtained from our numerical model. Chapter 6 introduces the 4DV AR method of data assim- ilation as well as the particular parametrizations of the physical drivers used in our implementation. Then in chapter 7 we show validation experiments for the 4DV AR approach. Chapter 8 contains additional details about both mathematical and pro- grammatic aspects of our implementation. Chapter 9 summarizes and concludes the thesis. 5 Chapter 2 Plasma physics of the ionosphere The Earth’s ionosphere is a layer of the atmosphere at altitude between 50 km to 2,000 km above the Earth’s surface. Although electrically neutral gases still dominate this region of atmosphere, the presence of free electrons and ions is the main focus of our interest in ionosphere. In particular, the interaction between radio waves and free electrons in the ionosphere leads to reflection of radio waves and delay in the wave propagation. As a result, the understanding of ionospheric physics plays a crucial role in the development of wireless communication systems. The fundamental quantities for ionospheric dynamics are the densities, temperature and velocity of neutral gases, ion species and free electrons [6], [19]. The interaction between the ions and free electrons with geomagnetic and electric fields is also an important aspect of the dynamics in ionosphere. The most important external driv- ing force for ionosphere is the solar radiation that ionizes the neutral gases. Thus, a complete model of the dynamics of ionosphere would include fluid dynamic equations for neutral and ion species, thermal dynamic equations for heat transfer between dif- ferent gas species, Maxwell equations for magnetic and electrical fields and chemical kinetic modeling the ionization and recombination processes. In general the following simplifying approximations are made: The relatively low number of ionized particles compared to the density of neutral gases makes the impact of ion species on the neutral gases negligible. As a result, most models of ionosphere treat neutral gases as independent of ionized 6 gases. That is, although density, temperature and velocity of the neutral gases play important roles in the fluid dynamic equations of the ion gases, the reverse interaction is absent in the equations for neutral gases. The changes of electric and magnetic fields due to the motion of free electrons and ions is sufficiently small compared to the changes of the magnetic field caused by solar activities. It is therefore more effective to consider both magnetic and electric fields as driving forces. The thermal balance of ionized gases is a direct function of solar angle and geomagnetic activity. Therefore, the temperatures of ionized gases are assumed to be known. Using these approximations, the main physical laws governing ionospheric dynam- ics are the law of balance of mass, also referred as continuity equation and law of bal- ance of mechanical momentum. These laws are written for each ionized gas separately. As a result, the balance of mass is equivalent to the balance of number density. For an ion speciesi2fO + ;H + ;He + ;N + 2 ;O + 2 ;NO + ;N + g, a state variablen i in our model is the number densities of the ion species. The physics equations in our model are written as follows. 2.1 Basic equations of an ionosphere model 2.1.1 The law of balance of mass The law of balance of mass states that the rate of change of the total number of particles of an ion species in a given volume V is equal to the sum of the net inward flux of 7 particles across the boundary ofV and the production, minus recombination of the ion particles inside of the volumeV . This is written as d dt Z V n i (t)dV = Z @V n i ~ v i ~ dS + Z V P i dV Z V R i dV; (2.1) or using the divergence theorem, @ t n i =r (n i ~ v i ) +P i R i (2.2) where, V is the volume element and @V is its surface; n i (t) is the density of ion i at timet; ~ v i is the velocity of the ion flux through the surface; ~ is the unit outward normal vector to the surface@V ;P i is the rate of production of the ion due to photo- ionization. The rate of ionization is proportional to the intensity of solar Extreme Ultra Violet (EUV) radiation, density of the neutral gas that is being ionized (i.e.,O forO + , He forHe + ,O 2 forO + 2 N 2 forN + 2 ). At a given altitude, the effective EUV radiation is reduced due to the absorption by neutral gases above. The optical depth at an altitude r for a given wavelength of radiation is given by (;r) = Z 1 r X j j ()n j (r 0 )dr 0 ; where j () is the photoabsorption rate of thej-th gas. The ionization rate is given by P i (r) = X () i ()n i e (;r) ; where the summation is over all relevant wavelength of solar radiation. 8 The termR i in (2.1) is the chemical recombination rate of the ioni and is given by R i = X j6=i ij n j ii n i (2.3) The first term in this equation represents the chemical production for the ioni and the second term represents the chemical loss for ioni. ij ’s are the chemical recombination coefficients given by an empirical model. Table 2.1.1 shows the chemical reactions modeled by our implementation. In general, a first order chemical reaction represented by the expression A +B!C +D is translated by the Rate Law into mathematical equations for the number densities of the species in the form d dt n C = n A n B ; d dt n A = n A n B ; where is the reaction rate. Thus, fori6= j, ij = P r r n r , where the sum is over all the reactions where thei-th ion was produced, andn r is the density of the neutral species in the reaction. Similarily, ii = P r r n r , where the sum is over reactions where thei-th ion was lost, andn r is the density of the neutral species or the electron density in the reaction. Remark 2.1.1. It is important to note here, that because of this, ii ’s depend on the electron density (and therefore on the densities of all the ions). But the coefficients for the chemical production, ij , wherei6=j do not depend on the ion densities. 9 Reaction Rate (m 3 s 1 ) 1) H + +O!O + +H 2:22T 0:5 i 2) He + +N 2 !N + 2 +He 3:5e 16 3) He + +N 2 !N + +N +He 8:5e 16 4) He + +O 2 !O + +O +He 8:0e 16 5) He + +O 2 !O + 2 +He 2:0e 16 6) N + +O 2 !NO + +O 2:0e 16 7) N + +O 2 !O + 2 +N 4:0e 16 8) N + +O!O + +N 1:0e 18 9) N + +NO!NO + +O 2:0e 17 10) O + +H!H + +O 2:5e 17 T 0:5 n 11) O + +N 2 !NO + +N 1:533e 18 5:920e 19 T i 300 + 8:600e 20 ( T i 300 ) 2 ; for 300KT i 1700K 2:730e 18 1:155e 18 T i 300 + 1:483e 19 ( T i 300 ) 2 ; forT i > 1700K 12) O + +O 2 !O + 2 +O 2:82e 17 7:74e 18 T i 300 + 1:073e 18 ( T i 300 ) 2 5:17e 20 ( T i 300 ) 3 + 9:65e 22 ( T i 300 ) 4 13) O + +NO!NO + +O 1:0e 18 14) O + +NO!NO + +O 1:4e 16 ( T i 300 ) 0:44 15) N + 2 +O 2 !O + 2 +N 2 5:0e 17 ( T i 300 ) 0:5 16) N + 2 +O 2 !NO + +NO 1:0e 20 17) N + 2 +NO!NO + +N 2 3:3e 16 18) O + 2 +N!NO + +O 1:2e 16 19) O + 2 +N!N + +O 2 2:5e 16 20) O + 2 +NO!NO + +O 2 4:4e 16 21) O + 2 +N 2 !NO + +NO 5:0e 22 22) H + +e!H 4:43e 18 ( Te 300 ) 0:7 23) He + +e!H 4:43e 18 ( Te 300 ) 0:7 24) N + +e!H 4:43e 18 ( Te 300 ) 0:7 25) O + +e!H 4:43e 18 ( Te 300 ) 0:7 26) O + 2 +e!O +O 1:8e 13 ( Te 300 ) 0:39 27) NO + +e!N +O 1:6e 13 ( Te 300 ) 0:55 28) N + 2 +e!N +N 4:2e 13 ( Te 300 ) 0:85 Table 2.1: Rates for the chemical reactions modeled in GAIM++ The rates and formulae are taken from SAMI2 ([17]). Here T i , T e and T n are the temperatures of the ions, electrons and the neutrals respectively. For an example, in reactions 1,4 and 8 in table 2.1.1, the ionO + is produced, and in reactions 10-14 and 25 the ionO + is lost. The chemical reaction termR O + thus has the form R O + = 1 n O n H + + 4 n O 2 n He + + 8 n O n N + ( 10 n H + 11 n N 2 + 12 n O 2 + 13 n NO + 14 n NO + 25 n e )n O +: 10 2.1.2 The law of balance of momentum The second important fluid dynamic law is the law of balance of momentum. In a given volumeV , the total mechanical momentum of the particles of thei-th ion species is given by Z V n i ~ v i dV: The law of balance of momentum states that the change of the total momentum is equal to the total forces acting on the ion particles in V . These forces include gas pressure, Earth’s gravitational force, forces created by collisions with neutral particles and other ion species and, forces of electromagnetic fields. One important simplify- ing assumption, the quasi-stationary assumption, is often made in the modeling of ionospheric dynamic. This assumption states that the fluid reaches a steady state in momentum nearly instantaneously, that is for almost any instancet, the rate of change of momentum is equal to zero. Under the quasi-stationary condition, the law of bal- ance of momentum becomes equation of balance of forces acting on the ion particles. This is given by: r(n i kT i ) +n i M i ~ g +c i n i ~ E +~ v i ~ B ~ W i = 0 (2.4) where the first termr(n i kT i ) represent the negative of pressure gradient, the second termn i M i ~ g represents the gravitational force, the third termc i n i ~ E +~ v i ~ B repre- sents the force of electromagnetic field and the last term ~ W i represents the ion’s force due to collisions with other ions and neutral particles. ~ W i is given by ~ W i =n i M i in (~ v i ~ u) + X j n i M i ij (~ v i ~ v j ); (2.5) 11 whereT i and ~ v i are the temperature and the velocity of thei-th ion, respectively. The summation in the above equation is over all ion speciesj. Other vector fields included in the above equations are: the gravitation field ~ g, the neutral wind velocity ~ u, the electric field ~ E and the geomagnetic field ~ B. The constants in these equations are: k is the Boltzmann constant; M i is the mass of the ioni; c i is the charge of ioni; in is the collision frequency of the ioni with neutrals; ij is the collision frequency of the ioni with the ionj. It is often assumed that the ion velocity perpendicular to the local geomagnetic field is primarily due to convective forces. It is possible to obtain a separate model for these components. In particular, in the low latitude equatorial region, the electrical field is primarily in the eastward longitudinal direction [1]. The ion drift velocity in ~ E ~ B direction is in the latitude plane. As a result, we write~ v i =v i;EB ~ ? +v i k ~ k , where ~ ? and ~ k are unit vectors in the geomagnetic longitudinal plane perpendicular and parallel to the local geomagnetic line, respectively. The component perpendicular to the geomagnetic field, v i;EB , is specified by an external ~ E ~ B velocity model. The component parallel to the geomagnetic field, v i;k , is derived from (2.4) and the Momentum Equation for electrons. The mass of a free electron is so small, that the momentum equation for free elec- trons becomes r(n e kT e ) +c e n e ~ E +~ v e ~ B = 0; (2.6) 12 where n e , T e , c e and ~ v e represent the number density, the temperature, the electric charge, and the velocity of free electrons, respectively. In our model, we assume elec- tric neutrality of the macroscopic medium. That is n e = X i n i : Taking the projection of (2.4) and (2.6) in the direction of ~ k , and using the fact that (~ v i ~ B) ~ k = (~ v e ~ B) ~ k 0, we obtain: r k (n i kT i ) +n i M i g k +c i n i ~ E ~ k ~ W i ~ k = 0; (2.7) and r k (n e kT e )en e ~ E ~ k = 0; (2.8) whereg k =~ g ~ k and the operatorr k is defined by r k (f) =r(f) ~ k : Equation (2.8) leads to ~ E ~ k = 1 en e r k (n e kT e ) = k e r k (T e ) k en e r k (n e )T e : (2.9) Thus by substituting ~ E ~ k into (2.7) we have r k (n i kT i ) +n i M i g k c i n i k e r k (T e )c i n i k en e r k (n e )T e ~ W i ~ k = 0: (2.10) 13 Letu k =~ u ~ k be the parallel component of the neutral wind velocity, then (2.5) becomes ~ W i ~ k =n i M i in v i k u k + X j n i M i ij v i k v j k (2.11) Combining (2.11) and (2.10) together we get r k (n i kT i ) +n i M i g k c i n i k e r k (T e )c i n i k en e r k (n e )T e =n i M i in v i k u k + X j n i M i ij v i k v j k (2.12) Our aim is to solve for the ion velocity parallel to the geomagnetic field linev i k . Equation (2.12) corresponds to a system of linear equations forv i k . This can be made more explicit by dividing both sides of equation n i M i and regrouping the terms, we obtain in v i k + X j ij v i k v j k = k n i M i r k (n i )T i c i k en e M i r k (n e )T e k M i r k (T i )c i k eM i r k (T e ) + in u k +g k (2.13) Since all the ions in our model have only one electron missing, thenc i =e and we will omit c i e from now on. LetJ be the number of ions and letn be the column vector of the ion densities, n = (n i ) i=1;:::;J and letv k be the column vector of velocities of ions parallel to geomagnetic line. Then (2.13) can be rewritten in matrix form as Av k =B(n)r k (n) +f 14 where A = 0 B B B B B B B @ 11 12 1J 21 22 23 J1 JJ 1 C C C C C C C A ; (2.14) where ii = in + P j6=i ij , B(n) =k 0 B B B B B B B @ 11 Te neM 1 Te neM 1 Te neM 1 Te neM 2 22 Te neM 2 Te neM 2 Te neM J Te neM J Te neM J JJ 1 C C C C C C C A ; (2.15) where ii = T i n i M i + Te neM i and f = k M i r k (T i ) k M i r k (T e ) + in u k +g k i=1;:::;J (2.16) Since all the ij ’s and in ’s are positive, this makes A diagonally dominant, and therefore invertible. Thus v k =A 1 B(n)r k (n) +A 1 f =D(n)r k (n) +w; whereD(n) =A 1 B(n) andw =A 1 f. Finally, we have v i k = X j d ij (n)r k (n j ) +w i : (2.17) 15 We note that the elements d ij (n) depend on densities of all the ions species. The following proposition proves that all the elements ofD are negative: d ij < 0;8ij (2.18) Proposition 2.1.1. LetA = (a ij ) be a matrix, with a ij 0; 8j6=i; P j a ij 0; 8i: 9 > = > ; (2.19) andAx =y. Ify> 0, thenx> 0. Proof. Without loss of generality letx 1 = min i (x i ). For allj6= 1, sincea 1j 0, this impliesa 1j x j a 1j x 1 . Then we have 0<y 1 =a 11 x 1 + X j6=1 a 1j x j a 11 x 1 + X j6=1 a 1j x 1 =x 1 X j a 1j and because P j a 1j 0, we getx 1 > 0, and thereforex> 0. According to (2.15) the columns ofB are negative. It is also clear from (2.14), that A satisfies the condition (2.19) of the proposition. Therefore, thej-th column ofD, D j =A 1 (B j ) is negative, which proves (2.18). 16 2.2 Complete equation for ionospheric dynamics in dif- ferential form Combining equation (2.17) for the ion velocity with the continuity equation (2.2), we obtain the complete differential form of ionospheric dynamic equation: @ t n i =r (n i v i ? ~ ? )r (n i v i k ~ k ) +P i R i : (2.20) By using the equality r (~ u) =r()~ u +r~ u; we obtain @ t n i = r ? (n i v i ? )n i v i ? r ~ ? r k (n i v i k )n i v i k r ~ k +P i R i : or, by pluggin in (2.17), @ t n i = r ? (n i v i ? )r k ( n i X j d ij (n)r k (n j ) +w i !) (2.21) n i v i ? r ~ ? n i X j d ij (n)r k (n j ) +w i ! r ~ k +P i R i : By now we have completely described the mathematical equations (2.21) which is the basis of the work in this thesis. We intentionally presented the equations in a coor- dinate free form. This allows a variety of selections of coordinates for the numerical implementation of the model. We also note that (2.21) represents a system of degen- erated nonlinear parabolic partial differential equations. In fact, the second derivative 17 of the density functions only appears in the direction parallel to the geomagnetic field line. These characteristics of the system present major challenges for both theoretical and numerical studies of these equations. 2.3 Models of the Earth’s geomagnetic field The Earth’s magnetic field varies in both intensity and direction as a function of geo- graphic location and time. Factors contributing to this variation including secular vari- ation induced by interaction of the liquid outer core of the Earth and geomagnetic field, residual magnetization of lithosphere, ionospheric current, solar activity. However, the dominant component of the Earth’s magnetic field, referred to as core field, changes very little as function of time. In this thesis, we approximate the Earth’s magnetic field by the core field. An empirical model of the core field, the International Geomagnetic Reference Field [23], is represented by a spherical harmonic potential function of the form (r;;;t) =R E N X n=1 R E r n+1 n X m=0 Q m n (;t)S m n (cos()); (2.22) wherer,, andt are distance from the center of the Earth, colatitude, longitude, and time; R E is the reference radius of the Earth, (R E = 6371:2km); S m n is the Schmidt semi-normalized Legendre functions of degree n and order m, and Q m n (;t) = (g m n (t) cos(m) +h m n (t) sin(m)): The Schmidt semi-normalized Legendre functions are defined by S m n (x) = (1) m s 2(nm)! (n +m)! P m n (x); 18 whereP m n is the Legendre function P m n (x) = (1) m (1x 2 ) m=2 d m dx m P n (x); andP n is the Legendre polynomial defined by P n (x) = 1 2 n n! d n dx n (x 2 1) n : The scalar valued function is referred to as the empirical geomagnetic potential function. That is, the Earth’s magnetic field vector ~ B is given by the negative gradient of . A geomagnetic field line is a 3-dimensional curve with tangent at every point on the curve equal to ~ B. As a result, we must integrate the following non-linear system of differential equations to find a geomagnetic field line: d ds ~ x(s) = ~ B(~ x(s)): In our ionospheric model, only the direction of the geomagnetic field is used. A dipole model of the geomagnetic field lines is adequate for an initial development of our model. A dipole model of the Earth’s ionosphere uses a potential field of the form q(r;) = R 2 E cos() r 2 : (2.23) The above potential field can be viewed as a truncated and scaled version of (2.22) withn = 1 andm = 0. Equation (2.23) is often referred to as a Earth centered dipole model of the Earth’s geomagnetic potential model. In a more general dipole model, the distance to the center of the Earthr is replaced by the distance to a point ~ E 0 away from the the center of the Earth and the colatitude is replaced by the colatitude in 19 an alternative coordinate system referred to as geomagnetic coordinate system. In the geomagnetic coordinate system the origin of the coordinate system is ~ E 0 and the z- axis is a line passing through ~ E 0 and tilted from the Earth’s rotational axis. The shift in the center of the dipole and the tilt in the axis can lead to a more accurate empirical dipole model of the Earth’s magnetic field. This type of model is often referred to as an eccentric-tilted dipole model. In the remainder of this thesis, we shall use a dipole model to approximate the Earth’s magnetic field. In most cases, an eccentric-tilted dipole model is used. In these cases, the variablesr and refer to the magnetic radius and magnetic colatitude. For a given value ofq, the curve with constant value ofq is a plane curve on the magnetic longitudinal plane. Another family of curves given by p(r;) = r R E sin 2 () ; (2.24) has the property that at every point (r;) the tangent vectors of the two curves, q- constant andp-constant, are perpendicular to each other. This can be verified by eval- uating the gradient vectors of the two functions. In fact, sincer = @ @r ~ r + 1 r @ @ ~ , we have rq(r;) = 2R 2 E r 3 cos()~ r R 2 E r 3 sin()~ ; rp(r;) = 1 R E sin 2 () ~ r 2 R E sin 3 () cos()~ ; where~ r and~ are unit vectors in the radial and angular direction, respectively (see figure 2.1. Sincerq andrp are perpendicular, we find that the tangent direction of the p-constant curve at any given point is in the direction of the gradient vector of the potential field q. As a result, we call a p-constant curve a geomagnetic field 20 line because the tangent of this curve gives the direction of the magnetic field. Simi- larly, aq-constant curve is referred to as a geomagnetic potential line. The families of geomagnetic field lines and geomagnetic potential lines form a orthogonal curvilinear coordinate system in a geomagnetic longitudinal plane. In the following section, we shall rewrite the coordinate free form of the ionospheric model equation as a system of partial differential equations in this coordinate system. μ p p-constant line μ q q-constant line Figure 2.1:p andq lines 2.4 Partial differential equations for ionospheric dynamics In equation (2.21), the gradient operatorr and its projectionsr ? andr k are in coor- dinate free form. When a dipole model of the Earth’s geomagnetic field is used, the 21 directions ~ ? and ~ k are the tangent directions of q-constant and p-constant lines, respectively. As a result, for any given functionf, we have r ? f(p;q) = @ p f(p;q)krpk; (2.25) r k f(p;q) = @ q f(p;q)krqk; (2.26) where krq(r;)k 2 = 4R 2 E r 6 cos 2 () + R 4 E r 6 sin 2 (); krp(r;)k 2 = 1 R 2 E sin 4 () + 4 R 2 E sin 6 () cos 2 (): We denotekrpk andkrqk bya p anda q , respectively. In order to translate the coor- dinate free form of the divergence operator into partial derivative with respect to the variablesp andq, we consider a small area in geomagnetic meridional plane bounded by two magnetic field lines withp-coordinatesp 1 andp 2 and two magnetic potential lines withq-coordinatesq 1 andq 2 . For any vector field ~ f, the divergence theorem gives Z r ~ f(p;q)dA = Z @ ~ f(p;q)~ dS; where~ is the outward normal of the boundary@ . Equivalently, we have Z p 2 p 1 Z q 2 q 1 r ~ f(p;q)a 1 p (p;q)a 1 q (p;q)dpdq = Z p 2 p 1 a 1 p (p;q 2 ) ~ f(p;q 2 )~ k (p;q 2 )a 1 p (p;q 1 ) ~ f(p;q 1 )~ k (p;q 1 ) dp + Z q 2 q 1 a 1 q (p 2 ;q) ~ f(p 2 ;q)~ ? (p 2 ;q)a 1 q (p 1 ;q) ~ f(p 1 ;q)~ ? (p 1 ;q) dq: 22 By dividing both sides of the above equation by (p 2 p 1 )(q 2 q 1 ) and taking the limit asp 2 p 1 andq 2 q 1 tend toward zero, we obtain r ~ f(p;q) = a p r k (a 1 p ~ f~ k ) +a q r ? (a 1 q ~ f~ ? ) = a p a q @ q (a 1 p ~ f~ k ) +a q a p @ p (a 1 q ~ f~ ? ): Using these equalities and (2.17), equation (2.2) can be written in the form of a partial differential equation as follows: @ t n i = a q a p @ p (a 1 q n i v i ? ) (2.27) a q a p @ q ( a 1 p n i X j d ij (n)a q @ q n j +w i !) +P i R i : For the model equation to be well-posed, we must also impose boundary conditions. In general, we consider a domain (see figure 2.2) of the form =fpp max andr r 0 g =f(p;q);p2 (p min ;p max );q2 (q 0 (p;r 0 );q 0 (p;r 0 ))g, where q 0 (p;r 0 ) = R 2 E r 2 0 r 1 r 0 R E p : For a small enough value of p min , the domain covers the equatorial region of the meridional plane above radius r 0 . The boundary @ = has two sides. The part of boundary 1 =fr = r 0 g =f(p;q);q =q 0 (p;r 0 ) orq = q 0 (p;r 0 )g corre- sponds to the lowest altitude of the domain. The upper part of boundary 2 =fp = 23 p max g =f(p;q);p = p max ;q2 (q 0 ;q 0 )g separates ionosphere from plasma sphere. The following boundary conditions are imposed. P i =R i ; for (p;q)2 1 ; (2.28) v i ? = 0; for (p;q)2 2 : (2.29) Σ 1 Σ 2 Ω p max p min r 0 Figure 2.2:pq region In practice, 1 is a surface with a specified geographic altitude. This leads to an interval inq that is not symmetric with respect to zero. For simplicity of our analysis, we assume here that 1 is a sphere centered at the origin of the geomagnetic coordinate system. We also note that (2.29) is not really a boundary condition on the solution. However, this condition guarantees that there are no ions escaping or entering from 24 the upper part of the boundary. In addition to these boundary conditions, we assume that the initial state of model is given by n i (t 0 ;p;q) =n i;0 (p;q): (2.30) Equations (2.27)-(2.30) completely define the basic model of the Earth’s ionosphere we are using in this thesis. 25 Chapter 3 Well posedness of equations In this chapter, we consider the well-posedness of the partial differential equations (2.27) in a domain of the form =f(p;q);p2 (p min ;p max );q2 (q 0 (p;r 0 );q 0 (p;r 0 ))g with boundary condition (2.28)-(2.29) and initial condition (2.30). It is interesting to note that (2.27) is a degenerated nonlinear parabolic equation since the second deriva- tive inp is missing on the right hand side of the equation. We first consider the ion velocity v i? in the direction perpendicular to the geo- magnetic field line. Consider a particle at a position (p(t);q(t)) in the geomagnetic meridional plane moving at a velocityv i? at a timet, the velocity can be written as v i? (t;p;q) =a 1 p (p;q)v ip (t;p;q); where v ip (t;p;q) = p 0 (t) is the rate at which the coordinate p changes in time. The rate of change of theq coordinate is equal to zero because the particle is moving in the tangent direction of aq constant curve. When the velocityv ip (t;p;q) is independent of the value ofq, all ion particles along a magnetic field line (p-constant line) have an identical rate of change in thep-coordinate. As a consequence, all ion particles along a magnetic field line with value p 0 at time t 0 will be found along another magnetic field line p 1 at a time t 1 . This synchronized motion in the p-direction gives rise to 26 the concept of ion flux tube. One can imagine that the ion particles along a magnetic field line are contained in a invisible tube that can shrink or expand in time to another magnetic field line. There is no exchange of particles among the flux tubes. This is indeed not the real situation in ionosphere. The velocityv ip (t;p;q) can, in general, be different for different values ofq. The difference in the rate of change ofp-coordinate of the ion position along a magnetic field line is referred to as the shear velocity along a magnetic field line. In most ionospheric models, a shear-free velocity field which correspond toq independentv ip (t;p;q) is assumed. In our numerical implementation of our ionospheric model, we do not assume shear-free velocity. However, we make a shear-free velocity assumption to introduce a significant simplification of our model. Under the shear-free velocity assumption, equation (2.27) can be written as @ t n i = v ip (t;p)@ p (n i ) +b(t;p;q)n i (t;p;q) (3.1) a q a p @ q ( a 1 p n i X j d ij (n)a q @ q n j +w i !) +P i R i ; where b(t;p;q) =a p a q @ p a 1 p a 1 q v ip (t;p) : We assume thatv ip are identical for all ion species and we defined the following char- acteristic curves int-p plane d ds t(s) = 1; and d ds p(s) =v ip (s;p(s)): (3.2) so that@ t n i v ip @ p n i =@ s n i . Using this definition, we can write (3.1) as 27 @ s n i = a 2 q @ q n i X j d ij @ q n j ! cn i X j d ij @ q n j (3.3) +a q w i @ q n i +f i n i +P i R i ; where c = a p a q @ q (a 1 p a q ) f i = b +a q @ q w i +a q a p @ q (a 1 p )w i : Under the shear-free velocity assumption, we have reduced a partial differential equation in 2-dimensional spatial domain (3.1) to a family of partial differential equa- tions in 1-dimensional spatial domain (3.3). Moreover, (3.3) is a non-degenerate parabolic equation. It is interesting to note that many numerical models of the Earth’s ionospheric solve (3.3) for a large number of flux tubes. More precisely, these model select a finite number of values ofp at timet 0 . Then the equation of characteristic curves (3.2) are solved for these initial values ofp. Using these solutions, (3.3) is solved along each characteristic curve. It is not difficult to see that as the flux tubes move in time, new flux tubes need to be introduced to cover the desired spatial domain. Similarly, when a flux tube moves outside the domain of interest, it is eliminated. The house keeping of these flux tubes represents a complex part of the model. Equation (3.3) represents a coupled system of nonlinear parabolic equations. In general, the well-posedness of these initial boundary value problems is not guaranteed. In the next section, we can discuss the well-posedness of the solution first for the single ion case and later for the case of multiple ion species. 28 3.1 Existence and uniqueness of the solution - Single ion case It is not uncommon that a model of the Earth ionosphere only describes density vari- ation of the predominant ion species,O + . In this case, the collisions betweenO + and other ion species are neglected. Such a model is referred to as a single ion ionospheric model. In a single ion ionospheric model, the computation of ion velocity parallel to geo- magnetic field line is greatly simplified. In fact, from (2.13) gives v ik = k(T i +T e ) in n i M i r k n i +w i : In the above equality, we used the fact that in a single ion model, we haven i = n e . In the remainder of this section, we omit the subscripti and assume all the quantities involved are forO + . The equation for ion density along a flux tube (3.3) becomes @ s n = a 2 q @ q k(T +T e ) M @ q n + c k(T +T e ) M +a q w @ q n (3.4) +fn +PR: It is easy to see that (3.4) is a linear parabolic partial differential equation with variable coefficient of the form @ s n =a 2 @ 2 q n +a 1 @ q n +a 0 n +h; (3.5) for allt> 0 andq2 (q 0 (p(s);r 0 );q 0 (p(s);r 0 )). Since a flux tube is moving in time, the range of values ofq of the flux tube at a specified altitude also changes as a function 29 of time. The condition at the lower boundary of a flux tube is given by the Dirichlet boundary condition n(s;q 0 (p(s);r 0 )) = n b (s;p(s);q 0 (p(s);r 0 )); n(s;q 0 (p(s);r 0 )) = n b (s;p(s);q 0 (p(s);r 0 )); where n b (s;p(s);q 0 (p(s);r 0 )) is a given function defined on a spherical surface at a distance r 0 from the Earth’s center. The above boundary condition is derived from chemical equilibrium condition at the radiusr 0 . The chemical equilibrium condition (2.28) assumes that the production and recombination rate at radiusr 0 are equal. This gives the ion density at the surface with radiusr 0 as a function of the production rate and recombination rate. The coefficientsa k in (3.5) are given by a 2 = a 2 q k(T +T e ) M ; a 1 = a 2 q k M @ q (T +T e ) + ck M (T +T e ) +a q w; a 0 = b +a q @ q w +a q a p @ q (a 1 p )w: The nonhomogeneous termh is given by h =P: We note that (3.5) is a linear parabolic partial differential equation defined on a time varying interval (q 0 (p(s);r 0 );q 0 (p(s);r 0 )). We can normalize the interval to 30 (1; 1) by replacingq withx =q=q 0 (p(s);r 0 ). That is, we define a function ^ n(s;x) = n(s;xq 0 (p(s);r 0 )). Then we have, @ s ^ n(s;x) = @ s n(s;xq 0 (p(s);r 0 )) +@ q n(s;xq 0 (p(s);r 0 ))@ p q 0 (p(s);r 0 )p 0 (s); @ x ^ n(s;x) = @ q n(s;xq 0 (p(s);r 0 ))q 0 (p(s);r 0 ); @ 2 x ^ n(s;x) = @ 2 q n(s;xq 0 (p(s);r 0 )) (q 0 (p(s);r 0 )) 2 : These equations lead to @ s ^ n(s;x) = @ s n(s;xq 0 (p(s);r 0 )) + 1 q 0 (p(s);r 0 ) @ x n(s;xq 0 (p(s);r 0 ))@ p q 0 (p(s);r 0 )p 0 (s): Therefore, (3.5) can be rewritten as @ s ^ n(s;x) = ^ a 2 (s;x)@ 2 x ^ n(s;x) + ^ a 1 (s;x)@ x ^ n(s;x) + ^ a 0 (s;x) ^ n(s;x) + ^ h(s;x); (3.6) where ^ a 2 (s;x) = a 2 (s;xq 0 (p(s);r 0 )) (q 0 (p(s);r 0 )) 2 ; ^ a 1 (s;x) = a 1 (s;xq 0 (p(s);r 0 ))@ p q 0 (p(s);r 0 )p 0 (s) (q 0 (p(s);r 0 )) ; ^ a 0 (s;x) = a 0 (s;xq 0 (p(s);r 0 )); ^ h(s;x) = h(s;xq 0 (p(s);r 0 )): 31 The boundary condition for ^ n is given by ^ n(s;1) = n b (s;p(s);q 0 (p(s);r 0 )); ^ n(s; 1) = n b (s;p(s);q 0 (p(s);r 0 )): If we define ~ n(s;x) = ^ n(s;x) 1x 2 ^ n(s;1) 1 +x 2 ^ n(s; 1); we have @ s ~ n(s;x) = ^ a 2 (s;x)@ 2 x ~ n(s;x) + ^ a 1 (s;x)@ x ~ n(s;x) + ^ a 0 (s;x)~ n(s;x) + ~ h(s;x); (3.7) where ~ h(s;x) = ^ h(s;x)+ 1x 2 d ds ^ n(s;1)+ 1 +x 2 d ds ^ n(s; 1)+ ^ a 1 (s;x) 2 (^ n(s;1)^ n(s; 1)): The boundary condition for ~ n is given by ~ n(s;1) = ~ n(s; 1) = 0: The scaling in q direction requires that q 0 (p(s);r 0 )6= 0 for all s. We make the following assumption. Assumption 3.1.1. TheEB drift velocityv p is a continuous function int andp and v p (t;p) equals to zero for allpp min andpp max . Moreover,q 0 (p max ;r 0 )q min > 0. 32 Lemma 3.1.1. Assumption 3.1.1 guarantees that if a flux tube has a non-zero length at a timet 0 , then for all timett 0 , the interval ofq values for the flux tube has a length larger than minfq 0 (p(t 0 );r 0 );q 0 (p max ;r 0 )g. Proof: The equation of characteristic (3.2) shows that d dt p(t) =v i;p (t;p(t)): Sincev i (t;p) = 0 forpp max andpp max , the interval (p max ;p max ) is invariant for the above dynamical system. As a result, for allp(t 0 )2 (p max ;p max ),q 0 (p(t);r 0 ) q min . Ifp(t 0 ) p max , since these flux tube do not move in time, q 0 (p(t);r 0 ) is also constant. The above result guarantees that the coefficients in (3.6) are bounded and smooth in time if the driving forces functions that are involved in these coefficients are smooth. In general, the following regularity conditions are satisfied by the driving forces of the ionosphere. Assumption 3.1.2. The driving forces such as the temperatures of ions and electrons, neutral wind velocity, neutral density and the collision frequencies are sufficiently smooth that the coefficients ^ a 2 , ^ a 1 , ^ a 0 and ~ h are inC 1 ((0;T ) (1; 1). Moreover, there exists a constant ^ a 2;0 > 0 such that ^ a 2 (t;x) ^ a 2;0 for allt> 0 andx2 (1; 1). In this case, using Theorem 2 in page 233 in [18], equation (3.6) has a unique solu- tion for any initial condition inC 4 (1; 1). The well-posedness condition of (3.6) can certainly be obtained with less restrictive conditions. However, this thesis is focused on the numerical approximation of (3.6), the approach of using finite difference approxi- mation to show the existence of the solution is also useful for the proof of stability of our numerical approximation. 33 3.2 Existence and uniqueness of the solution - multiple ion case In the previous section, we make use of the fact that the single ion case of the model for ionosphere reduces the model equation to a linear parabolic partial differential equation. In the case of a multiple ion model of the ionosphere, we must work with a complex nonlinearly coupled system of parabolic partial differential equation. In fact, along a single flux tube, a multiple ion species equation has the form: @ t ~ n(t;x) =A 1 (t;x;~ n)@ 2 x ~ n(t;x) +A 2 (t;x;~ n)@ 1 x ~ n(t;x) +A 0 (t;x;~ n)~ n(t;x) + ~ h(t;x): (3.8) Even though (3.8) is semilinear, the coupling of the second order term makes the well- posedness of (3.8) challenging to establish. In the subsequent Chapter of this thesis, we can establish the uniform boundedness of the solution of (3.8). 34 Chapter 4 Formulation of numerical approximation In the previous chapter, we have shown that the system of partial differential equations governing the dynamics of the ion densities in the Earth’s ionosphere can be trans- formed into a family of nonlinear parabolic partial differential equations in one space dimension. In the case of single ion equation, we have shown that the equation is well-posed. However, it is extremely difficult to develop any analytical form of the exact solution in either single or multiple ion case. In this chapter, we shall formulate a numerical scheme to approximate the solutions of the partial differential equations established in Chapter 2 and 3. The method of characteristics that transforms the partial differential equations in three space dimensions to a family of partial differential equations along moving flux tubes has also led to a numerical approach that approximate the solutions of the family of parabolic differential equations. This approach is widely used in the ionospheric modeling community (see [2], [3], [5] and [4]). It is often referred to as the Lagrangian frame approach because the 1-dimensional modeling grid follows the motion of the flux tube. The primary goal for our development of a new model for ionosphere is to use the resulting model for data assimilation. The tracking of moving flux tubes not only introduces challenging house-keeping issues, it also makes it difficult to propagate the model uncertainty matrix in time. As a result, we have chosen an approach in which 35 the modeling grid is fixed relative to the geomagnetic coordinate system. This type of approach is commonly referred to as the Euler frame approach. The basic technique we are using to derive our numerical method is called the finite volume approach [10], [7], [11]. This approach is motivated by the integral form of the model equations. After subdividing the modeling domain into small volume elements (also called ”volume-pixels” or voxels), we approximate the densities of ions inside a voxel by a constant function. Then, the basic laws of balance of mass and mechanical momentum are written for each voxel. This leads to an approximation of the partial differential equations similar to the method of finite difference. In Section 4.1, we describe the subdivision of the the modeling domain. In Section 4.2, we formulate the finite volume approximation of the model equations. 4.1 Space discretization In chapter 2 we have described the (p;q;) coordinate system, where the p- and q- coordinates are good approximations for the Earth’s Geo-Magnetic Field and Geo- Magnetic Potential lines. Recall that q = R 2 E sin r 2 (4.1) p = r R E cos 2 (4.2) We note that the model equations correspond to a degenerated system of parabolic equations. In particular, there is no diffusion in the direction perpendicular to the geo- magnetic field line. In order to take advantage of this separation of physical properties in these orthogonal directions, we choose to have voxels that are rectangular shaped boxes in (p;q;)-coordinate system. 36 As a result, our space discretization consists of choosing appropriate grids for each coordinate, such that it covers the region of interest; has sufficient spatial resolution in areas where the ionosphere has most structure, that is in the low altitudes, as well as low to mid latitudes, and the size of the elements changes smoothly when transitioning from higher to lower resolution areas. Figure 4.1 shows that a simple uniform grid in (p;q) does not satisfy these require- ments. We can see from (4.2) that the altitude resolution at the magnetic equator (where = 0) is a function ofp only. Therefore, the selection of the grid inp is determined by our selection of the resolution in altitude. More precisely, the following strategy is followed. First p min is chosen so that the intersection of the p-line with the mag- netic equator corresponds to the point of the region with the lowest altitude. On the other hand p max is chosen to be the smallest p such that the entire region of interest is included in the area delimitted by the p-line (figure 4.2). In order to control the grid and make it denser in the lower altitudes that higher altitudes, another p-value, p midlat , is selected to be thep-value at the point of intersection of a pre-specified lat- itude (called the transition latitude, usually about 30 degrees) and the highest altitude boundary of the region of interest. A uniform grid is created betweenp midlat andp max using a specified p step size. To create a denser grid between p min and p midlat , a smooth function is chosen to transition from p to a smaller step size in the lower altitudes. This algorithm of choosing the grid inp, particularly the choice ofp min and p max , shows that our model region will always be larger than the region of interest. 37 4 6 8 x 10 6 −8 −6 −4 −2 0 2 4 6 8 x 10 6 Figure 4.1: Uniform grid inp andq We note that the spacing in q determines the latitude resolution in low and mid latitudes, as well as, the altitude resolution in high latitudes. As a result, we need to balance the needs for creating appropriate latitudinal grid in low and mid-latitude region and altitude grid in high latitude region when we select a grid inq for a global model. However, since the model presented in this thesis is primarily for low and mid latitude, our selection of grid in q is determined by the need for appropriate latitude resolution. A low resolution grid used in GAIM++ is shown in figure 4.3. The decision to use a rectangular grid in (p;q;)-coordinate system has its own downside as well. In 38 p max p midlat θ max p min r min r max θ trans Figure 4.2: Minimum, maximum and midlat values ofp general, the horizontal feature size in ionosphere is much larger than the vertical fea- ture size. In particular, vertical resolution in the F-region of the ionosphere (around 200-400km altitude) needs to be lower than 10km to accurately capture the character- istic of the ionosphere. On the other hand, horizontal resolution of several hundred kilometer is considered acceptable. The use of a rectangular grid creates a grid that couples the vertical resolution with the horizontal resolution. In this thesis, we work with a rectangular grid in (p;q;)-coordinate system only. However, an alternative space discretization strategy can be explored in the future. 39 4 6 8 10 x 10 6 −8 −6 −4 −2 0 2 4 6 8 x 10 6 MX (km) MZ (km) Figure 4.3: A low resolution GAIM++ grid Letp 1 < p 2 < < p max andq 0 < q 1 <q max are the grids thus chosen. And let 1 < 2 << max is a uniform grid in magnetic longitude, covering the region of interest. The volume element defined by this grid is given by V =f(p;q;) :p n p<p n+1 ;q m q<q m+1 ; k < k+1 g (4.3) 40 In addition to the volume elements, the other geometric elements also important for the definition of the approximation of the model equations are the facets of the voxels. In fact, each voxel has 6 facets defined by @V p = f(p;q;) :p =p n ;q m q<q m+1 ; k < k+1 g @V p + = f(p;q;) :p =p n+1 ;q m q<q m+1 ; k < k+1 g @V q = f(p;q;) :p n p<p n+1 ;q =q m ; k < k+1 g @V q + = f(p;q;) :p n p<p n+1 ;q =q m+1 ; k < k+1 g @V = f(p;q;) :p n p<p n+1 ;q m q<q m+1 ; = k g @V + = f(p;q;) :p n p<p n+1 ;q m q<q m+1 ; = k+1 g As we mentioned in section 3.1 the lower boundary of our region is handled using the chemical equilibrium condition. Therefore we need to identify the lower boundary elements of our discretized region, on which this condition will need to be set in our numerical equation. Definition 4.1.1. A voxel V of our grid is called inactive if either one of it’s neighbor- ing voxels is below our region. All the other voxels are called active. Figure 4.4 shows the inactive voxels on our GAIM grid which form the lower boundary for our region. 4.2 Numerical form of the equations The basic approach for the approximation of the integral form of the balance of mass equation is the finite volume method. The elementary volumeV is a voxel defined in 41 4 6 8 10 x 10 6 −8 −6 −4 −2 0 2 4 6 8 x 10 6 MX (km) MZ (km) Figure 4.4: Inactive voxels (in blue) defining the lower boundary of the GAIM++ grid the previous section. Letk be the discrete time index, and t k =t k+1 t k be thek-th time step. The left hand side of the equation (2.1) is approximated (at timet k ) by: d dt Z V n i (t)dV! n k+1 i n k i t where! is the volume of the voxel. 42 Given a voxel from the grid described in section 4.1, the surface integral terms of the right hand side of the equation (2.1) are approximated by the product of the surface area and the ion flow rate across the faces of the voxel as follows: Z @V n i ~ v i ~ dS X f=p ;q ; f f v if n k if ; (4.4) where f is the face index (f =p ;q ; ); f is the “scaled“ area of facef of the voxel; f is 1 for lower surfaces, and1 for upper surfaces; n k if is the density of the ioni on the center of the facef at timet k ; v if = ~ v i ~ f is the velocity of the ioni at the center of the facef. Notice, thatv if has a instead of a time index. This is because it will be treated differently forq-faces. In fact, as we have discussed in Chapter 3, the ion velocity in the direction perpendicular to the magnetic field lines is specified by an empirical model. As a result, the ions are carried by this convective force across the boundary. In this case, the flow rate is determined by the product of ion velocity and the up-stream ion density. Forp- and -faces however the velocity will be treated with the explicit scheme in time, and the ion density across the surface is computed using ion densities inside of voxel at timek. The two families of surfaces parallel to the magnetic field lines are those with constant values and those with constantp values. The normal directions of these surfaces are perpendicular to magnetic field lines. In low and mid-latitude, the ion drift velocity is negligible in the zonal direction which is normal to the constant surfaces. As a result, the empirical model forEB velocity sets the zonal velocity to zero. However, in the polar region, the zonal drift is the predominate component ofEB velocity. In this thesis, our discussion is focused primarily on mid and low latitude region of 43 the ionosphere. Thus, we will now omit the terms corresponding to ion flow across constant surfaces from our equation. In contrast to surfaces parallel to the magnetic field lines, the normal direction of aq-surface is parallel to the magnetic field line. As we have discussed in Chapter 2, the ion motion is driven by a combination of diffusion and convection. This is char- acterized by the fact that the ion velocity across aq-surface is determined, in part, by the difference of ion density on two sides of the surface. On the other hand, the grav- itational force, the collisions with neutral particles and the effect of the temperature gradient are among the convective forces that carry ions acrossq-surface. As a result, a combination of explicit scheme and implicit scheme in time is used to evaluate the flow rate acrossq-surfaces. In particular, we use primarily implicit scheme to compute the diffusion terms and explicit scheme for the convective terms. The use of an implicit scheme is motivated by the stability of the scheme. Since in the multi-ion model the diffusion term is nonlinear, we use an explicit scheme to evaluate the coefficients of the diffusion term while use an implicit scheme to evaluate the ion density gradient in the normal direction of theq-surface. More precisely, to discretizev iq , we notice that our choice of the orthogonal (p;q) coordinate system, makesv q =v k . So we use (2.17), to get: v iq =d ii r q (n iq ) + X j6=i d ij r q (n jq ) +w i (4.5) Here the first term will be treated with the implicit scheme, while the rest will be treated with the explicit scheme: v iq = d k+1 ii r q (n k+1 iq ) + X j6=i d k ij r q (n k jq ) +w k i = d k+1 ii r q (n k+1 iq ) + ~ v k iq (4.6) 44 where ~ v k iq = X j6=i d k ij r q (n k jq ) +w k i (4.7) The approximation of the volumetric integrals in (2.1) is given by Z V P i dVp k i ! (4.8) wherep k i is the production rate in the voxel for ioni at timet k . And, finally, using (2.3), the chemical reaction term is approximated by Z V R i dV k ii n k+1 i X j6=i k ij n k j ! ! (4.9) where n k i is the density of the ioni in the center of voxel at timet k ; k ij is the chemical recombination coefficient at timet k . We note that in (4.9) the chemical production term is treated with the explicit scheme, while the loss is treated with the implicit scheme, even though the loss rate k ii uses the explicit scheme. Combining all of the above, we obtain ! n k+1 i n k i t = X f=q f f v if n k if + X f=p f f v k if n k if +p k i ! k ii n k+1 i ! + X j6=i k ij n k j ! 45 Substitutev if discussed previously into (4.10), multiply both sides by t ! , to get n k+1 i n k i = t X q=q q ! q v iq n k iq + t X f=p f ! f v k if n k if + tp k i t k ii n k+1 i + t X j k ij n k j = t X q=q q ! q d k+1 ii r q (n k+1 iq ) + ~ v k iq n k iq + t X f=p f ! f v k if n k if + tp k i t k ii n k+1 i + t X j k ij n k j (4.10) Taking all the implicit terms to the left-hand side of the equation, and the explicit terms to the right-hand side, we obtain n k+1 i + t k ii n k+1 i t X q=q q ! q d k+1 ii n k iq r q (n k+1 iq ) =n k i + t X f=p f ! f v k if n k if + t X q=q q ! q ~ v k iq n k iq + t X j6=i k ij n k j + tp k i (4.11) The last two quantities for which we need to specify the discretization are the gra- dient of the density and the densities on the faces. For the gradient, let the " superscript indicate the voxel on theq-upper side of the face and # the index of theq-lower side of the face. Then r q (n q +) = n " n d " r q (n q ) = nn # d # 46 where n " and n # are the densities in the neighboring voxels on the q-upper and q- lower sides, andd " andd # are the distances between the voxel and its neighbors on the q-upper andq-lower sides. For the density on the facen if in the right hand side of the equation, it is taken to be the density of the voxel where particles move out from. So we look at the direction of the velocity vector on that face (v if or ~ v iq forq-faces). If this velocity is directed inward to the voxel, thenn if is set to be the density of the neighboring voxel. Otherwise, if the velocity points outward, then the density on the face is set to be the velocity in the voxel. For the density on the face in the left hand side of the equation (and every other place they are used), the density is set to be the average of the densities of the two neighboring voxels. In these specifications, the case in which the voxel is a boundary voxel, where the needed neighboring voxel does not exist, several strategies can be considered. For example, for the gradient we consider it to be 0 (i.e. the density in the non-existing voxel is the same as the density in the voxel). In equations (4.11) we still have to choose t appropriately. Once t is chosen, these equations can be written in the following matrix form An k+1 =Bn k +P (4.12) wheren k is the vector of all densities for all ions in all the voxels, i.e. densities of each ion in the first voxel, followed by densities of each ion in the second voxel, etc. Thus, if we haveN voxels, andJ ions, thenn k andF have dimensionsNJ1,A andB have dimensionsNJNJ. Also, bothA andB are clearly sparse matrices, withA having at most 3 non-zero elements per row, andB having at most 6 +J non-zero elements per row. If we also number the voxels in the order ofq-neighbors (i.e. q-neighboring 47 voxels succeed each other in the linear order), thenA will also be a banded (in this case tri-diagonal) matrix. An even more important property of A is that it is diagonally dominant, more precisely we have a ij 0; 8j6=i; P j a ij > 1; 8i: 9 > = > ; (4.13) Indeed, we note that in equation (4.11),r q (n q ) is multiplied by q . Taking into account (2.18), n k+1 in the sum on the left hand side, whose coefficient adds to the diagonal term ofA, enters the equation with a positive sign from both termsq andq + . At the same time the densities in the neighboring voxels,n " andn # , whose coefficients form the off diagonal terms of A, both enter it with negative signs. Thus the first property is clearly verified. Also, by the same observation, the contribution of the terms in that sum to the off diagonal terms is exactly the same as their contribution to the diagonal term, except with the opposite sign. Therefore P j a ij = 1 + ii , and since ii > 0, this proves the second property. The diagonal dominance ofA guarantees that it is invertible. Moreover, thanks to (4.13), by applying Proposition 1 toA, we have that ifx > 0, thenA 1 x > 0. This means that if we make sure that for every timek,Bn k +P > 0, then we will always get a positive density state at timek + 1 as a solution to (4.12). This sets a constraint on how to choose t in (4.11): t has to be small enough such thatBn k +P> 0. This will give an upper boundary on t. We will continue on the choice of t in section 4.3. Remark 4.2.1. Equation equation (4.12) holds only on the active voxels (see definition 4.1.1). This means that the rows ofB andP corresponding to the inactive voxels are set 48 to 0. Then on each step, after the equation is solved and the next density state is found, the corresponding elements of that state are set by using the chemical equilibrium formula given by (2.28). 4.3 Time discretization and Stability analysis Earth’s ionosphere is extremely dynamic, the ion densities change very rapidly in time. Thus, even if our equations allowed for large t’s we would still restrict to some values where we can consider the state of the ionosphere fairly static. In GAIM++ implementation we typically set an upper limit of 12 minutes on t. As mentioned above, to insure that our density state stays positive, we would need to choose t so that the right hand side of (4.12) is positive. (4.11) shows the form of that right hand side: B k n k +P k = (1 + t k V k + t k C k )n k + t k P k (4.14) whereV k represents the terms with ion velocities (2nd and 3rd terms in (4.11)), C k represents the term with chemical reactions (4th term in (4.11)), andP k is the vector of solar production. From this form it is very easy to constraint t k as to insure the positivity of the sum. We now proceed to the stability analysis of our equation, which will give us the final recipe on choosing t k . To guarantee the stability of our system, we need to show that the state vectorn k does not grow exponentially withk. More precisely we will prove that given an initial state vectorn 0 > 0,k( Q N k=1 A 1 k B k )n 0 k does not grow exponentially withN. In what follows, we will ommit the indexk when the statement applies to allk’s and when there is no ambiguity. 49 We start withA and prove thatA 1 is a contracting transformation. That is Proposition 4.3.1. Given (4.13), all the eigenvalues ofA are> 1. Proof. Let be an eigenvalue, with x the corresponding eigenvector, so that Ax = x. Without loss of generality let x 1 = max i (x i ) > 0 (otherwise we would have consideredx, with an appropriate permutation). Since a 1j 0, this implies a 1j x j a 1j x 1 for allj6= 1. Then x 1 =a 11 x 1 + X j6=1 a 1j x j a 11 x 1 + X j6=1 a 1j x 1 =x 1 X j a 1j >x 1 and becausex 1 > 0, we get> 1. We now would like to estimatekBnk 1 , wheren> 0. We start with the second term in (4.14) by noting that tVn represents the movement of ions from some voxels to others within the time t. It also models the “falling off” of some ions out of the lower boundary of our region. The fact that we don’t have production of ions coming in from the lower boundary, shows that the total mass of the ions represented byknk 1 = P l n l cannot increased due to ion movements. That is, the sum of the elements ofn + tVn is not larger thanknk 1 . So, if in addition we also require that 1 + tV 0 (4.15) then we getkn + tVnk 1 knk 1 . 50 The 3rd term in (4.14),Cn, models the production of ions due to chemical recom- bination reactions. From (4.11), we can see thatC is a block-diagonal matrix, with each block having the form 0 B B B B B B B @ 0 12 1J 21 0 2J J1 J2 0 1 C C C C C C C A whereJ is the number of ions. As we described in section 2.1.1, for i6= j, ij ’s depend on the temperatures of the ions and the temperatures and densities of the neutral species (see remark 2.1.1). So, since they do not depend on the ion densities, then there is a constant , such that 0 < ij < for all i6= j = 1;:::;J and for all voxels. Then it’s clear that kCnk 1 (J 1) knk 1 . Combining all of the above, we get kBnk 1 = kn + tVn + tCnk 1 kn + tVnk 1 + tkCnk 1 knk 1 + t(J 1) knk 1 = (1 + (J 1) t)knk 1 . Together with proposition 4.3.1 above, this gives k(A 1 B)nk 1 kBnk 1 (1 +ct)knk 1 wherec is a constant. Therefore, after theN-th step N Y k=1 A 1 k B k ! n 0 1 N Y k=1 (1 +ct k ) ! kn 0 k 1 (1 +ct) 1 t kn 0 k 1 kn 0 k 1 where t = max k=1;:::;N (t k ). This proves the stability of our system. 51 The only requirement necessary to achieve this was set by (4.15). Intuitively, (4.15) means, that after a round of shuffling the ions, no voxel will end up empty (or in deficit of ions), no matter the state we started off from. The inequality (4.15), gives a trivial method on choosing t on each step, that would guarantee both the positivity of our solution and the stability of our system. Our choice of the implicit-explict scheme and the algorithm for choosing t, allows us to preserve the positivity of our density state while keeping the numerical algorithm stable. However, the effect of this on t is that it could become rather small especially in the multiple ion case when used with high resolution grids. This presents practical computational challenges. There are several possible remedies to this issue, such as parallelizing the algorithm so it can be run in a clustered environment; using a nested-grid approach to restrict the high resolution grids to small regions, etc. Changes to the numerical scheme, where more terms can be treated with the implicit scheme (and thus relax the condition resulting from (4.15)) are also being considered. 52 Chapter 5 Validation of the forward model An extremely important part of the development of a model of the Earth’s ionosphere is the validation of the model. Since the model predicted ion densities are used as a prior estimate for ion density for data assimilation, it is critical that the model generated ionosphere density distribution has both the known characteristics of the ionospheric ion density and the expected sensitivity to ionospheric driving forces. In this chapter, we present an extensive array of numerical results from our model aimed at estab- lishing the validity of our model. In Section 5.1, we focus our attention at the main characteristics of the model output. These characteristics include vertically integrated density or Vertical Total Electron Content (VTEC) in both magnitude and spatial dis- tribution, vertical profiles of ion species, as well as, cross sections of the ion density in meridional planes. In Section 5.2, we present the sensitivity of the model density with respect to driving forces. This is particularly important for our work in the estimation of ionospheric drivers. Naturally, a more direct comparison between measured ion density and model pre- diction is desirable for model validation. However, direct measurements are extremely sparse in both spatial and temporal domain. On the other hand, in the course of the use of our model in data assimilation, we can indirectly validate our forward model. 53 5.1 Characteristics of the model ionosphere The most widely-known characteristic of the ionosphere is the vertical profiles of the ion densities. In particular, the vertical distribution of the 3 most important ion species O + ;H + ;He + . It is known that the O + is the most dominant ion in altitude lower than 1,000km with a peak density at around 200-300km altitude. The peak ofO + ion coincides with F-region peak. However, at altitude above 1,000kmH + becomes the dominant species. In addition to the general shape of the vertical profile, the magnitude of the peak density is also an important characteristic of the ionosphere. In the Figures 5.1 and 5.2, vertical profiles of all three ion species are shown. The main characteristics of these profiles are highly consistent with published results of similar profiles, such as the ones shown in figure 5.3 from [4]. 54 10 6 10 8 10 10 10 12 0 500 1000 1500 Altitude (km) LT: 06:00, Lat: 5, Lon: 0 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 LT: 18:00, Lat: 5, Lon: 180 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 Altitude (km) LT: 06:00, Lat: 30, Lon: 0 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 LT: 18:00, Lat: 30, Lon: 180 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 Density (1/m 3 ) Altitude (km) LT: 06:00, Lat: 55, Lon: 0 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 Density (1/m 3 ) LT: 18:00, Lat: 55, Lon: 180 O + H + He + Figure 5.1: GAIM++ ion density profiles for O + , H + , He + at dawn and dusk for different latitude regions as of Sep 21, 2001. Note that 2001 is a high solar activity year 55 10 6 10 8 10 10 10 12 0 500 1000 1500 Altitude (km) LT: 06:00, Lat: 5, Lon: 0 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 LT: 18:00, Lat: 5, Lon: 180 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 Altitude (km) LT: 06:00, Lat: 30, Lon: 0 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 LT: 18:00, Lat: 30, Lon: 180 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 Density (1/m 3 ) Altitude (km) LT: 06:00, Lat: 55, Lon: 0 O + H + He + 10 6 10 8 10 10 10 12 0 500 1000 1500 Density (1/m 3 ) LT: 18:00, Lat: 55, Lon: 180 O + H + He + Figure 5.2: GAIM++ ion density profiles for O + , H + , He + at dawn and dusk for different latitude regions as of Sep 21, 2006. Note that 2006 is a low solar activity year 56 Figure 5.3: Ion density profiles forO + (—),H + (– – –),He + (- - -) from [4]. Note that here the ion densities are given incm 3 , whereas in the plots for our model they are inm 3 Another widely known characteristic of the ionosphere is the phenomenon called equatorial anomaly. This refers to the observation that the vertical total electron con- tent has lower value at or near the equator than the VTEC value in mid and low latitude regions after 2PM local time. This phenomenon is caused by the process of ion foun- tain. In fact, as the solar radiation intensity increases the ion density near the equator through photo-ionization along magnetic field lines with low value ofp, the diffusion along these field lines is limited. At the same time, the ions are pushed by EB forces upward into high altitude magnetic field lines with higher value of p. At the higher altitude, ion diffusion velocity is much higher than at low altitude. The ions diffuse toward lower altitude part of these flux tube with higher value ofp. The result- ing distribution of ions exhibits two separate high density regions on either side of the equator. In figures 5.4, 5.5, 5.6 ion densities in a meridional plane at different local 57 times are shown. One can clearly see the on set of photo-ionization following sun rise. The bulk part of high ion density region is over the equatorial region at this phase. Later during the day, we see a rise of high density region in altitude. This followed by a downward diffusion along magnetic field lines creating two separate high density regions. It is worth noting that the ionospheric production rate which is a function of neutral density and solar radiation intensity does not have this spatial split. The model recreation of the equatorial anomaly indicates that our forward model correctly captures the convective and diffusive process of the ionosphere. 58 500 1000 1500 LT 09:00 (45 o ) Altitude (km) Density 0.5 1 1.5 2 2.5 3 3.5 x 10 12 LT 12:00 (90 o ) 500 1000 1500 LT 15:00 (135 o ) Altitude (km) LT 18:00 (180 o ) −50 0 50 500 1000 1500 LT 21:00 (225 o ) Latitude (deg.) Altitude (km) −50 0 50 LT 24:00 (270 o ) Latitude (deg.) Figure 5.4: GAIM++O + ion density longitude slices 59 500 1000 1500 LT 09:00 (45 o ) Altitude (km) Density 1 2 3 4 5 6 7 x 10 10 LT 12:00 (90 o ) 500 1000 1500 LT 15:00 (135 o ) Altitude (km) LT 18:00 (180 o ) −50 0 50 500 1000 1500 LT 21:00 (225 o ) Latitude (deg.) Altitude (km) −50 0 50 LT 24:00 (270 o ) Latitude (deg.) Figure 5.5: GAIM++H + ion density longitude slices 60 500 1000 1500 LT 09:00 (45 o ) Altitude (km) Density 1 2 3 4 5 6 x 10 9 LT 12:00 (90 o ) 500 1000 1500 LT 15:00 (135 o ) Altitude (km) LT 18:00 (180 o ) −50 0 50 500 1000 1500 LT 21:00 (225 o ) Latitude (deg.) Altitude (km) −50 0 50 LT 24:00 (270 o ) Latitude (deg.) Figure 5.6: GAIM++He + ion density longitude slices 61 A more direct validation of the forward model is through comparison of VTEC maps generated by our model with purely data driven reconstruction of VTEC using ground based GPS measurements. In fact, at the JPL group, global VTEC maps (Global Ionospheric Map, GIM) are generated hourly using slant TEC data derived from measurements of signal delays between ground receivers and GPS satellites. GIM has been extensivly calibrated and validated against other independent measure- ments of the Earth’s ionosphere. In figures 5.7 and 5.8 our forward model generated VTEC map for solar min and solar max conditions are shown. The general character- istics of these VTEC maps are consistent with the purely data driven GIM shown in figures 5.9 and 5.10. 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 20 40 60 80 100 120 140 Figure 5.7: GAIM++ Vertical TEC map as of 2001/09/21 06:00 UT, solar maximum year 62 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 6 12 18 24 30 36 42 Figure 5.8: GAIM++ Vertical TEC map as of 2006/09/21 06:00 UT, solar minimum year 90°S 60°S 30°S 0° 30°N 60°N 90°N 0° 60°E 120°E 180° 120°W 60°W 0° Longitude (deg) Latitude (deg) 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 120.0 130.0 140.0 Figure 5.9: GIM Vertical TEC map as of 2001/09/21 06:00 UT, solar maximum year 63 90°S 60°S 30°S 0° 30°N 60°N 90°N 0° 60°E 120°E 180° 120°W 60°W 0° Longitude (deg) Latitude (deg) 0.0 10.0 20.0 30.0 40.0 Figure 5.10: GIM Vertical TEC map as of 2006/09/21 06:00 UT, solar minimum year 64 5.2 The sensitivity of the model to some physics drivers One of the primary motivations for the development of the new version of the iono- spheric forward model is to use this model in the estimation of ionospheric driving forces. It is therefore extremely important for us to examine the sensitivity of the model output to the perturbation of these driving forces. In particular, since our most abundant observation is ground based slant TEC measurements, we use the VTEC maps generated by our new model to examine the sensitivity of our model prediction to the changes of the driving force. 5.2.1 ExB Drift One of most often used characteristics of the ionosphere in the estimation of driving force is the distance between the two high density regions of the equatorial anomaly. This distance is correlated to the strength of EB drift velocity. As the EB velocity increases, ions are pushed to increasingly higher flux tubes. This leads to a larger distance between the two high VTEC regions around the equator. In figures 5.12 and 5.13, we can see that our model outputs correctly reproduced this dependency. 65 0 5 10 15 20 −60 −40 −20 0 20 40 60 80 100 Local time (hrs) ExB velocity (m/sec) original perturbed Figure 5.11:EB drift pre-reversal enhancement perturbation 10 8 10 9 10 10 10 11 10 12 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Density (1/m 3 ) Altitude (km) O + original O + ExB perturbed Figure 5.12: The sensitivity ofO + to the perturbation of theEB drift pre-reversal enhancement. Location: 5 N, 60 E, 06:00 UT 66 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 20 40 60 80 100 120 140 (a) No perturbation 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 20 40 60 80 100 120 140 (b) EB perturbed Figure 5.13: GAIM Vertical TEC map, with EB drift pre-reversal enhancement perturbation 67 5.2.2 Neutral Wind Another driving force affecting the equatorial anomaly is the presence of the neutral wind. In particular, the direction of neutral wind along a magnetic field line. As in any one-dimensional convection-diffusion equation, the presence of convection skews the symmetric process of diffusion. For an example, a strong north-ward neutral wind gives preference to northward diffusion relative to southward diffusion of ions along a flux tube. This leads to higher density and VTEC north of the equator relative to the high density region south of the equator. Figures 5.15 and 5.16 show the effect on the equatorial anomaly for different strength of neutral wind. We can see the enhancement of non-symmetry as neutral wind velocity increases. Local Time [hrs] Latitude [deg] 2 4 6 8 10 12 14 16 18 20 22 −25 −20 −15 −10 −5 0 5 10 15 20 25 −10 0 10 20 30 40 50 60 70 80 90 100 Figure 5.14: Wind perturbation g 68 10 8 10 9 10 10 10 11 10 12 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Density (1/m 3 ) Altitude (km) O + original O + Wind perturbed Figure 5.15: The sensitivity ofO + to the perturbation of neutral wind. Location: 10 N, 75 E, 06:00 UT 69 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 20 40 60 80 100 120 140 (a) No perturbation 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 20 40 60 80 100 120 140 (b) Neutral wind perturbed Figure 5.16: GAIM Vertical TEC map, with Wind perturbed in the Northern hemi- sphere 70 5.2.3 Solar Production The most dramatic change in ionosphere is, of course, the change of production rate. However, ion production rate is not in general considered as an independent driving force because it is a combined effect of solar radiation and neutral density. As we can see from the discussion in Chapter 2, neutral density affects many aspects of iono- spheric dynamics. On one hand, the neutral density determines the optical depth that reduces the intensity of solar radiation. On the other hand, neutral density also affects chemical reaction rates. In our modeling approach, we allow direct estimation of ion production rate without separating the contributions from solar radiation and neutral density changes. This choice allows us to directly change the ion density by changing production rate. In Figure 5.17 and 5.2.3, we can see that the changes in production lead to the desirable effects in the change of VTEC. 10 8 10 9 10 10 10 11 10 12 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Density (1/m 3 ) Altitude (km) O + original O + Production perturbed Figure 5.17: The sensitivity ofO + to the perturbation of solar production. Location: 5 N, 30 W, 06:00 UT 71 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 28 56 84 112 140 168 196 (a) No perturbation 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 28 56 84 112 140 168 196 (b) Solar Production perturbed Figure 5.18: GAIM Vertical TEC map, production increased 1.4 times 72 Chapter 6 The 4DV AR Process The variational data assimilation technique considers the partial differential equations governing a dynamical system as a constraint when minimizing the difference between observations collected over a temporal and spatial domain. Instead of adjusting the estimation of the system state at a timet based only on observations made at or near timet and a model generated prior estimation of the system state for timet, variational data assimilation also uses observation considered over a period of time to adjust model parameters. The inclusion of spatial and temporal dimension of the observations leads to the term 4-Dimensional Variational Data Assimilation or 4DV AR [29]. In most common cases, the model parameter of interest is the initial condition of the system. In these cases, the 4DV AR technique chooses an initial state of the system so that the forecast of the system state provides the most consistent prediction to observations collected over a period of time. The practical implementation of 4DV AR is often asso- ciated with the use of the adjoint method to obtain efficient evaluation of the gradient vector of the optimization cost functional. However, in our application of 4DV AR, our parameters are the driving forces of the ionosphere. This is due in part to the strong diurnal cycle of the Earth’s ionosphere. The largest source of uncertainty in the dynamics of the Earth’s ionosphere is the driving forces. 73 6.1 The 4DV AR equations For the simplicity of the presentation, in this section we will restrict ourselves to the single ion case. With that in mind, we can see from (3.1), that the equation for our single ion model has the form @ t n =a@ 2 q n +b@ q n +c@ p n (6.1) where a represents the diffusion rate, b represents the convection rate, and c = c(p) represents the EB drift. Notice that since c only depends on p, we only need to know it on the equator. Now, assume thata,b andc are parametrized with a parameter vector, i.e. given the initial staten(t 0 ) and given, our forward model can compute the staten(t;). Given a set of observations of this model, the objective of 4DV AR is to find the parameter, that best “fits” these observations. As we saw in Chapter 4, the equation (6.1) is time-discretized based on the physics of the model. This type of discretization leads to a partitiont 0 ;t 1 ;:::;t n of the full time interval, where the segments [t k1 ;t k ] are not necessarily equal in length. Let’s say we get a set of observations distributed in the full time interval. Since the time of each observation is unrelated to our time discretization, each of the intervals [t k1 ;t k ] can receive a different number of observations (including 0). We will assume that all those are concentrated at timet k and we will call the vector of these concen- trated observationsy k (which could be of size 0). Consequently, the [t k1 ;t k ] interval should not be too large, since otherwise, the observations which have fallen into that interval would have a large variance and the state at timet k alone would not be able to represent them meaningfully. 74 In Chapter 4 we described the process of discretizing the equation (6.1) into the matrix equation A k n(t k+1 ) =B k n(t k ) +P k (6.2) LetH k be the observation matrix at timet k corresponding to observationy k . We form the following cost function J() = n X k=1 ky k H k n(t k ;)k 2 +k ^ k 2 (6.3) where ^ is an a priori stabilization parameter which we don’t want our solution to diverge too much from, and is the corresponding a priori weight. Thus, more formally, the goal of 4DV AR, is to find the that minimizes the cost functionJ(). To solve this minimization problem, we need to estimate r J = 0 B B B B @ @J @ 1 . . . @J @q 1 C C C C A The adjoint method consists of an efficient way of estimatingr J. Taking the partial derivatives in (6.3), we have @J @ j = n X k=1 2 y k H k n k T H k @n k @ j + 2 ( ^ ) T e j (6.4) wheree j = (0;:::; 1;:::; 0) T andn k =n(t k ;). 75 The only unknowns in these equations are @n k @ j . So we take the partial derivative of (6.2) with respect to j @ @ j A k n k+1 = @ @ j B k n k +P k @A k @ j n k+1 +A k @n k+1 @ j = @B k @ j n k +B k @n k @ j + @P k @ j Or by regrouping A k @n k+1 @ j =B k @n k @ j +F kj (6.5) where F kj = @A k @ j n k+1 + @B k @ j n k + @P k @ j (6.6) Notice, that since the dependence of A k , B k and P k on is known analytically, then their partial derivatives with respect to j and therefore F kj can be computed analytically. The equation (6.5) is called the sensitivity equation. Since the initial state n 0 = n(t 0 ) does not depend on, we also have @n 0 @ j = 0 (6.7) Even though the equations (6.5) and (6.7) would allow calculating @J @ j , it would be very inefficient, since we would need to solve the forward equation (6.5)M times (whereM is the number of parameters). We define the innovation vectorg k by g T k =2(y k H k n k ) T H k (6.8) 76 so, the first term in the right hand side of the equation (6.4) becomes @ ^ J @ j = n X k=1 2 y k H k n k T H k @n k @ j = n X k=1 g T k @n k @ j (6.9) We now choosev k so that 8 > < > : A T k1 v k1 B T k v k = g k ; fork = 1;:::;n v n = 0 (6.10) Equation (6.10)-(6.10) is called the adjoint equation. Notice that sinceA k is invert- ible, given theg k ’s, all thev k ’s can be computed recursively, starting fromn and mov- ing backwards. Using this, (6.9) becomes @ ^ J @ j = n X k=1 g T k @n k @ j = n X k=1 A T k1 v k1 B T k v k T @n k @ j = n X k=1 v T k1 A k1 @n k @ j n X k=1 v T k B k @n k @ j = n1 X k=0 v T k A k @n k+1 @ j n X k=1 v T k B k @n k @ j (k!k + 1) = n1 X k=0 v T k A k @n k+1 @ j n1 X k=1 v T k B k @n k @ j (sincev n = 0) = n1 X k=0 v T k A k @n k+1 @ j n1 X k=0 v T k B k T @n k @ j (since @n 0 @ j = 0) = n1 X k=0 v T k A k @n k+1 @ j B k @n k @ j = n1 X k=0 v T k F kj (from (6:5)) 77 So finally @J @ j = n1 X k=0 v T k F kj + 2 ( ^ ) T e j (6.11) This allows as to skip solving the sensitivity equation (6.5) for each element of the parameter vector, and only solve the adjoint equation (6.10) once (per estimation iteration). Afterr J is computed, a standard iterative algorithm is used to minimize J and find . For example, GAIM++ uses the well known LBFGS [9] algorithm to find the minimizing. 6.2 4DV AR in GAIM++ In GAIM++, the 4DV AR estimation technique is used for adjusting three physics drivers: the solar production,EB drift and the neutral wind. We now examine in more details the parametrization of each of these drivers, and derive the formulae for the terms needed to computeF kj in (6.6). 6.2.1 Parametrizations for the drivers Gaussian smoothing Let W = (W 1 ;:::;W N ) be some predefined set of fixed points inR 3 . Let be a function defined on these points. We can extend toR 3 by interpolating it with Gaussian smoothing: (X) = N X j=1 (W j )e d 2 (X;W j ) ;8X2R 3 whered is some distance function. 78 Thus, ifV = (V 1 ;:::;V K ) is another set of fixed points inR 3 , then (V l ) = N X j=1 (W j )e d 2 (V l ;W j ) ;8j = 1;:::;K This can be written in matrix form, as (V ) =T(W ) (6.12) where T = (T lj ) = e d 2 (V l ;W j ) is a matrix of dimension KN, and defines a smooth mapping fromW toV . Gaussian perturbation for solar production in local time Suppose now, that we have a predefined Earth-fixed spherical grid, and suppose W 1 ;:::;W N represent the centers of the voxels in that grid. Let = ( 1 ;::: N ) = ((W 1 );:::;(W N )). Then, assumingV 1 ;:::;V K are the centers of the voxels from the PQL grid described in section 4.1, the smooth mapping in (6.12) allows us to project onto ourPQL grid. This mapping allows us to parametrize the production driver, which is defined on our largePQL grid, by which is defined on a relatively coarse spherical grid. However, since ourPQL grid is Earth-fixed, we need to adjust the mapping to take into account the fact that the production is local time dependent. Thus, for a given timet, we’ll use the mapping (V ) =Z(t)T(W ) =M(t)(W ) (6.13) whereZ(t) = (z lj (t)) is the rotation matrix for timet, andM(t) = (m lj (t)) =Z(t)T . 79 Let ^ p(t;V l ) be the a priori production for one of the ions in the voxelV l at timet (see equation (4.8)). Then we set p(t;V l ;) = ^ p(t;V l ) 1 + X j m lj (t) j ! 2 Note: the reason of using parametrization of type(1+) 2 is to keep the production positive for any parameter. It is clear then, that @p(t;V l ;) @ j = ^ p(t;V l )2m lj (t) 1 + X j m lj (t) j ! It is easy to see from (4.11), that in our matrix equation (4.12) A and B don’t depend on the productionP , andP k = t k P k . Therefore, from (6.6) F kj (V l ) = @P k (V l ) @ j = t k @p k (V l ;) @ j = t k ^ p k (V l )2m lj (t k ) 1 + X j m lj (t k ) j ! Gaussian perturbation for neutral wind in local time The parametrization of the neutral wind is very similar to that of the production, since it’s using a similar mapping from a coarse spherical grid to the PQL. There are 2 differences however. First, since our forward model uses the neutral wind on the q faces, the pointsV 1 ;:::;V l in this case will represent the centers of theq faces. Second, since for the wind we don’t have a positivity constraint, the actual parametrization is a bit simpler, namely, if ^ u(t;V l ) is the a priori neutral wind velocity on the faceV l at timet, then we set u(t;V l ;) = ^ u(t;V l ) + X j m lj (t) j 80 and therefore @u(t;V l ;) @ j =m lj (6.14) From (4.11), we see thatA andP don’t depend on the neutral velocity u, and it entersB from the term with ~ v q . ~ v q in turn depends onu via (4.7) and (2.16), and one can see from these equations, that the elements ofB have a linear dependency onu. Therefore @B @ j can be calculated using (6.14) above. We will omit the calculation of the actual formula for it, as it is extremely long and is of no interest here. Perturbation for ~ E ~ B drift in local time In our forward model the ~ E ~ B drift velocity is based on an empirical model which gives the velocities on Erath’s magnetic equator, and then extends it to any point (r;;)2R 3 by setting v p (r;;) =v EB (r;;) =v eq () cos 3 p 1 + 3 sin 2 where v eq is the equatorial vertical drift velocity, and depends only on the longitude (or, equivalently, the local time). We use this representation to parametrize the ~ E ~ B drift velocity for the 4DV AR estimation. Namely, we conaider the following perturbation to the equatorial drift. If ^ v eq is the a priori equatorial drift, then we set v eq () = ^ v eq () + N X j=1 j ' j () 81 where' 1 ;:::;' N are some periodic “basis” functions on [0; 2]. In our implementa- tion we use a set of periodic cubic spline functions as this basis (see Figure 6.1). With this parametrization, we have @v p (r;;) @ j =' j () (6.15) Similar to the neutral wind, (4.11) shows thatA andP don’t depend onv p , and the dependency ofB’s cells is linear. Therefore @B @ j can be calculated using (6.15) above. Longitude (deg) 0 50 100 150 200 250 300 350 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 6.1: Cubic spline basis functions forEB drift perturbation In Chapter 8 we give more details about the algorithm and the implementation of 4DV AR in GAIM++. 82 6.2.2 4DV AR difficulties and issues Relative scaling of the drivers The physics drivers that we are adjusting with 4DV AR have different impacts on the model, the sensitivity of the model to each of them is different. If all of the drivers, or at least two of them are being adjusted at the same time (normally all should be adjusted), this can present difficulties for the minimization algorithm, since the land- scape of the cost function will be very non-linear (it will have very steep valleys). The get rid of higher order of non-linearities in the cost function before calling the mini- mization routine, we need to scale the drivers relative to each others. GAIM++ has a configuration option for this type of scaling. In pratice the actual scaling factors are either determined experimentally, or can be computed and adjusted at runtime based on the cost function partials. Dependency on t We said in section 6 above that the dependency of A, B and P on the parameter is known analytically, and thus their partial derivatives can be calculated. In practice however, this can be difficult, if these matrices have a dependency on the parameter through the time step t that the forward solver takes, because t may not be easy differentiable with respect to the parameter. As an example, if we examine howB is composed in section 4.2, we see that it clearly depends on t. On the other hand the strategy of choosing t in section 4.3 implies that it depends on the ion velocities. These velocities in turn are functions of theEB drift and the neutral wind, which we perturb in our 4DV AR process. Moreover, this dependency is not “smooth”, since it involves the min function. Currently, GAIM++ ignores these dependency. This is actually not a problem as far as theEB drift is concerned, because the value of t 83 is dominated by thev k velocity, which does not depend onEB. It does however depend on the neutral wind velocity. Thus, in the case of adjusting the neutral wind, ignoring this dependency introduces some error in our calculation of the gradient of the cost function. Non-linearity of the equation In our forward equation, the left hand side matrix A depends on the current state, and therefore depends on any parameter. That however is never accounted for in our 4DV AR equations or implementations. Currently, GAIM++ does not handle these dependencies, and considers the partial ofA with respect to all three parameters to be 0. On the other hand, in the right hand side, theB matrix also depends on the state, because of the v k ’s dependency on the gradient of the density state through the col- lisions between different ions (see (4.7). So again, there is a dependency that is neglected by the 4DV AR equations, and by GAIM++ implementation. Note how- ever, that this dependency does not exist in the single ion case, since there are no ion collisions to take into account. It is a subject of future research, to either account for all these dependencies in our computations, or have an estimate of the error produced as result of neglecting each or any of them, and see whether it can impact the minimization process. However, as can be seen from the OSSE results below, even with all the errors introduced by these simplifications, our implementation still works remarkably well. 84 Chapter 7 Validation of the 4DV AR Method The validation of a data assimilation model is an extensive process of cross compar- ison with reference data and tuning of model parameters. Substantial performance improvements can be obtained through this process. An earlier version of GAIM has undergone several years of validation efforts in a collaborative effort between JPL and USC teams. The results of these efforts lead to performance that represents the state of the art in ionosphere data assimilation ([21]). The earlier efforts are primarily based on recursive estimation techniques such as Kalman filter. Limited experience was also gained for 4DV AR. In this Chapter, we present our validation of 4DV AR approach using the new model. However, since the predominant ion in the region of ionosphere of interest is O + , we run our new model in a single ion mode in our evaluation of 4DV AR model. The results presented in this Chapter represent the first phase of validation of a data assimilation model. In this phase of validation, synthetic observations generated using our forward model is used. 7.1 OSSE When developing data assimilation algorithms, Observing System Simulation Experi- ment, or OSSE is an important tool for studying many different aspects of the system. In our case, it can help to test the correctness and evaluate the effectiveness of the 85 4DV AR process described above. But more importantly it allows us to evaluate differ- ent ways of parametrizing the drivers, understand the effects of estimating more than one driver at a time, and explore the retrievability of the parameters depending on the type and the quality of the observation data used. The experiment consists of 4 steps a) Generate an artificial (but “physical”) ionosphere. This is done by choosing a set of parameters (called the “true” parameter) for one or more of the drivers from section 6.2.1, applying the corresponding perturbations, then running the forward model. The ion density states thus generated is called “the wheather” and is supposed to represent the “true” ionosphere. b) Simulate observations. A real GPS STEC observation dataset is used as a basis for simulating the observations corresponding to the artificial ionosphere above. To better approximate the real world situation, a random white noise may be added to the observations. c) Assimilate the synthesized observations. This consists of running the 4DV AR process starting off from a “wrong” density state (called “the climatology”) and using the simulated observations above to adjust the drivers by finding the parameter minimizing the cost function. d) Assess the impact of the synthesized observations. This consists of comparing the reconstructed density states to the “wheather”, and the minimizing parameter to the “true” parameter. Below, we will present several results from OSSE runs, where the three drivers are perturbed and estimated alone and together in some combinations. The following settings were used for all the experiments below 86 We have used a low resolution low latitude region between30 and 30 degrees as our geometry. The reason to use a low latitude region is that the meridonal EB drift is only predominant in that region. For step b) above, slant TEC observations were simulated. The time and geom- etry of the observation links were based on real GPS TEC observation data. A random white noise with standard deviation of 2:0 TEC Units was added to the simulated observations. The same parametrization used for simulating the observations was also used for the estimation step c) above. This means, that the true parameter could be retrieved exactly if there was no noise added to the data. Eventhough this is no longer true when noise is added, we should still be able to recover a good approximation of the true parameter. For step c) above, the ionosphere density state generated by running the forward model with no perturbation is taken to be the “climate”. When running the 4DV AR program, cycles of 4 hours were used with a 1 hour overlap. This results in 8 cycles for a 24 hour period. See section 8.6 for more detailed explanation of cycles and the iteration process. We will use several “metrics” to asses the success of each experiment. First, we will look at the decrease of the cost function during the entire iteration process. When considering the cost function, it is important to keep in mind that a new portion of the observation data is used for each consecutive cycle. Because of this, eventhough a cycle starts off using the last state from the previous cycle, it is possible for the cost function to jump up compared to its value from the previous cycle. As the cycle 87 progresses however, the cost function should decrease. In fact, only the last iteration in a cycle, where the cycle is considered converged, matters for our estimation: the parameter we get at the end of each cycle is considered to be the minimizing solution for our problem for that cycle. Second, we will compare the prefit and postfit residuals by looking at their his- tograms. In our case, prefit residuals are the differences between the TEC observations corresponding to the climate and to the wheather. Similarly, the postfit residuals are the differences between observations corresponding to the state after a converged cycle and the observations corresponding to the wheather. We expect the distributions for these residuals to go from a widely spread to a narrow one approaching the distribution of the noise we added to the data (in our case a normal distribution with mean 0 and standard deviation 2.0). Third, we will examine the behavior of the parameter itself during the iteration. We expect it to start from 0 (since our climate is generated with no perturbation) and to approach the true parameter as the iterative process progresses. In the case ofEB drift, instead of showing the perturbation parameter we will show the actual equatorial drift pattern as it gives a better intuition about the convergence. We start by 3 experiments, where a single driver is perturbed, then estimated alone. 7.1.1 Solar Production For this experiment, to generate the “wheather” state of the ionosphere, the solar pro- duction was perturbed by uniformly increasing it by 21%. Thus a single parameter was being estimated. In terms of the parameterization of the perturbation described in section 6.2.1, this corresponds to having one big volume element in the spherical 88 grid. The physical meaning of this parameter can be described as the underlying solar radiation intensity. 0 3 6 9 12 15 18 21 24 10 4 10 5 10 6 Cost Function Time (hrs) and Cycles Figure 7.1: Cost function decrease over the 8 cycles when estimating solar production −20 −10 0 10 20 30 40 50 60 70 80 0 0.05 0.1 0.15 0.2 0.25 TECU normalized count prefit, mean: 18.57, stddev: 16.91 cycle 8, mean: −0.06, stddev: 2.04 Figure 7.2: Prefit residuals (bottom) and postfit residuals after the end of the last cycle 89 0 5 10 15 20 25 30 Iterations in cycles Increase in production (%) climate wheather iteration end of cycle Figure 7.3: The convergence of the production parameter. All the iterations in all the cycles are shown. The end of iterations are marked in red. 7.1.2 ExB Drift For this experiment, the “wheather” is generated using a perturbation that enhances the pre-reversal of the EB equatorial drift. The perturbation uses 9 parameters, corresponding to the 9 cubic basis spline functions shown in figure 6.1. 90 0 3 6 9 12 15 18 21 24 10 4 10 5 10 6 Cost Function Time (hrs) and Cycles Figure 7.4: Cost function decrease over the 8 cycles when estimatingEB drift −70 −60 −50 −40 −30 −20 −10 0 10 20 30 0 0.05 0.1 0.15 0.2 TECU normalized count prefit, mean: −15.85, stddev: 13.76 cycle 8, mean: −0.03, stddev: 2.04 Figure 7.5: Prefit residuals (bottom) and postfit residuals after the end of the last cycle 91 0 3 6 9 12 15 18 21 24 −60 −40 −20 0 20 40 60 80 100 local time (hrs) ExB velocity (m/sec) climate weather estimation cycle 2 iter 12 (a) End of cycle 2 0 3 6 9 12 15 18 21 24 −60 −40 −20 0 20 40 60 80 100 local time (hrs) ExB velocity (m/sec) climate weather estimation cycle 8 iter 4 (b) End of cycle 8 Figure 7.6: TheEB pattern convergence 7.1.3 Neutral Wind As described in section 6.2.1, a coarse spherical grid is used to parametrize the neutral wind driver. In this experiment we used a grid of dimensions 1x6x24 which results in 92 144 parameters being estimated. The actual perturbation is shown in figure 5.14, and consists of increasing the wind velocity in the northern hemisphere between 18 and 22 hours local time. 0 3 6 9 12 15 18 21 24 10 4 10 5 Cost Function Time (hrs) and Cycles Figure 7.7: Cost function decrease over the 8 cycles when estimating neutral wind −40 −30 −20 −10 0 10 20 30 40 0.1 0.05 0 0.05 0.1 0.15 0.2 TECU normalized count prefit, mean: 1.66, stddev: 5.05 cycle 8, mean: −0.03, stddev: 2.03 Figure 7.8: Prefit residuals (bottom) and postfit residuals after the end of the last cycle 93 When examining figure 7.8, it is important to note, that unlike the other cases seen so far, the distribution of the prefit residuals is not spread as widely (as compared to the noise we added to the data), nor is its mean very far away from 0. However, as can be seen from the other figures in this section, even this type of challenging situation is handled successfully by the 4DV AR process, and the perturbation pattern is recovered pretty well (see figure 7.9(b)). 94 longitude latitude 0 50 100 150 200 250 300 −30 −25 −20 −15 −10 −5 0 5 10 15 20 −10 0 10 20 30 40 50 60 70 80 90 100 (a) End of cycle 2 longitude latitude 0 50 100 150 200 250 300 −30 −25 −20 −15 −10 −5 0 5 10 15 20 −10 0 10 20 30 40 50 60 70 80 90 100 (b) End of cycle 8 Figure 7.9: The neutral wind pattern convergence We now proceed with examining cases where a combination of two drivers is per- turbed and estimated together. 95 7.1.4 Solar Production and ExB Drift together In this experiment both the solar production and theEB drift are perturbed at the same time, each using the perturbation described in sections 7.1.1 and 7.1.2 respec- tively. 0 3 6 9 12 15 18 21 24 10 4 10 5 10 6 Cost Function Time (hrs) and Cycles Figure 7.10: Cost function decrease over the 8 cycles when estimating solar production andEB drift together −50 −40 −30 −20 −10 0 10 20 30 40 50 0 0.05 0.1 0.15 0.2 0.25 TECU normalized count prefit, mean: 7.16, stddev: 18.81 cycle 8, mean: −0.01, stddev: 2.04 Figure 7.11: Prefit residuals (bottom) and postfit residuals after the end of the last cycle 96 0 3 6 9 12 15 18 21 24 −60 −40 −20 0 20 40 60 80 100 local time (hrs) ExB velocity (m/sec) climate weather estimation cycle 2 iter 12 −5 0 5 10 15 20 25 Increase in production (%) (a) End of cycle 2 0 3 6 9 12 15 18 21 24 −60 −40 −20 0 20 40 60 80 100 local time (hrs) ExB velocity (m/sec) climate weather estimation cycle 8 iter 8 −5 0 5 10 15 20 25 Increase in production (%) (b) End of cycle 8 Figure 7.12: The production andEB drift pattern convergence 7.1.5 ExB Drift and Neutral Wind together In this experiment both theEB drift and the neutral wind are perturbed at the same time, each using the perturbation described in sections 7.1.2 and 7.1.3 respectively. 97 0 3 6 9 12 15 18 21 24 10 4 10 5 10 6 Cost Function Time (hrs) and Cycles Figure 7.13: Cost function decrease over the 8 cycles when estimatingEB drift and neutral wind together −50 −40 −30 −20 −10 0 10 20 30 0 0.05 0.1 0.15 0.2 TECU normalized count prefit, mean: −10.98, stddev: 11.97 cycle 8, mean: 0.11, stddev: 2.79 Figure 7.14: Prefit residuals (bottom) and postfit residuals after the end of the last cycle 98 0 3 6 9 12 15 18 21 24 −60 −40 −20 0 20 40 60 80 100 local time (hrs) ExB velocity (m/sec) climate weather estimation cycle 8 iter 35 Figure 7.15:EB pattern at the end of cycle 8 longitude latitude 0 50 100 150 200 250 300 −30 −25 −20 −15 −10 −5 0 5 10 15 20 −150 −100 −50 0 50 100 150 Figure 7.16: Neutral wind pattern at the end of cycle 8 99 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 20 40 60 80 100 120 140 (a) Climate 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 20 40 60 80 100 120 140 (b) Wheather 0° 60°E 120°E 180° 120°W 60°W 60°S 30°S 0° 30°N 60°N 60°S 30°S 0° 30°N 60°N 0 20 40 60 80 100 120 140 (c) End of cycle 8 Figure 7.17: Vertical TEC maps for theEB and neutral wind estimation 100 The cost function in figure 7.13 shows that the convergence is not as good as in previous cases. However the histograms in figure 7.14 tell us that the residuals were reduced to the noise level. By examining the patterns in figures 7.15 and 7.16 we can see that the experiment was not entirely successful. According to these figures, it looks like the process “prefers” to adjust the neutral wind rather thanEB drift, which results in the drift pattern not moving much from its climate position, and the wind pattern converging to some unphysical pattern. Figure 7.17 in turn shows that the process did converge “in the vertical TEC space”. All this suggests that by tuning the penalty coefficients in the cost function, it would be possible to make the process adjust the drift more. And indeed, our experiments show that by heavily penalizing the wind term we can get a good convergence for the drift, however in that case the wind pattern does not converge well. Note, that even if it was possible to find specific values of these settings as to make both of the drivers converge properly to their true values, those settings would be problem dependent, and would not necessarily be useful in general, nor especially for the real data runs (as opposed to assimilated data). Thus we don’t currently claim to have a definitive answer as to why the process does not converge to the expected result for this case. Besides the obvious possibility of issues in our implementation code, there are many other possible explanations, among which - we may have not found the correct set of settings for this type of runs. These include the values for penalties in the cost function (see (6.3)), the scaling of the parameters (see section 6.2.2), and many other settings that can be tuned in GAIM++ 101 - the signatures of the model sensitivity to theEB drift and to the neutral wind perturbation might be too similar (at least for the perturbations we’ve used), that makes retrieving them together very difficult. - the ground based TEC observations may not be appropriate for retrieving the drift and neutral wind together. - our parametrization of these drivers is not appropriate for an estimation of both of these drivers together. It is part of ongoing research to investigate this case, and find a strategy for this type of estimation. During the course of this project we have performed many other types of OSSE- like experiments, most of them successful. Below we give short descriptions of some of them, without much details, because of the the limitted scope of this work. Simulate observations without adding any noise. In this case we were able to retrieve the true parameters exactly. Perturb one driver for simulating the observations, but estimate two of them. Perturb two drivers for simulating the observations, but estimate only one of them. Perturb the f10.7 value for simulating the observations and estimate the solar production (with one parameter). Test many different knots configurations for the spline basis functions in the EB parametrization. This includes using one set for the perturbation and a different set for estimation. 102 Test many different settings for the penalty term when estimating the neutral wind. Test many different configurations for the cycles, with and without overlaps. 103 Chapter 8 Design and Implementation of GAIM++ 8.1 Basic concepts 8.1.1 Memory and other resource management GAIM++ uses a few techniques to avoid explicit memory/resource management and yet be efficient and exception-safe. This includes the use of uBLAS [31] libraries for vector and matrix representation, and the use of STL [28] for other containers (of course other components of STL are used as well). The concept of RAII (Resource Acquisition Is Initialization) is also used extensively. Whenever dynamic memory allocation is necessary, boost::shared ptr’s [12] are used instead of raw C++ pointers. These techniques, mainly because of the use of C++ templates, add some syntactic noise to the code, and if unfamiliar, a developer may have a hard time filtering out the “real” code, which contains the algorithm or logical flow. However, this approach also reduces the clutter and the errors that can occur when doing the ”classical” style explicit memory management (as an example, notice that there is not a single delete statement in GAIM++!). 104 8.1.2 Input/Output To persists or load data about the vectors and matrices used by different modules, GAIM++ uses some specialized file formats, such as MATLAB [22] or NetCDF [30]. This has the advantage of being able to easily use data that is produced by GAIM++ in many 3rd party applications (Matlab, octave [27], ...). The inverse of this is also useful, when data produced by outside tools, can be directly fed into GAIM++. To abstract out the actual format in which the data is stored, GAIM++ uses a small hierarchy of ”matrix file” classes (see figure 8.1) AbstractMatrixFile MatlabFile NetcdfFile MatrixFile Figure 8.1: MatrixFile class hierarchy in GAIM++ MatrixFile is a pure virtual class, it defines a minimal interface to be imple- mented by the subclasses. That interface has functions to load or save matrices repre- sented by raw pointers. AbstractMatrixClass uses these primitive routines, to provide many useful functions on how to load/save different types of special matrix or vector containers, such as sparse matrices, banded matrices, 3-dimensional arrays, etc. 105 Thus for a particular format (e.g. Matlab), the implementing class (MatlabFile) can just implement the trivial interface fromMatrixFile, and then it will automat- ically be able to load/save higher level constructs mentioned above. So, by using the AbstractMatrixFile, the code is file format agnostic, and matrix IO operations are extremely simple. This design also makes it easy to support a different file format when needed. GAIM++ also includes some utility classes (SparseMatrixFile and SparseMatrixFieldFile) to serialize sparse matrices into a binary format, which is specifically designed to allow efficiently appending new submatrices to an existing one, and to efficiently retrieving any submatrix of the full one. This is use- ful for serializing the observation matrices used in the 4DV AR program (see below), because the same observation file needs to be read multiple times during the process. 8.1.3 Loader/Generator apparatus On top of the I/O abstraction discussed above, GAIM++ provides another layer of abstraction for accessing the model related data. The idea is to releave the higher level numerical algorithms from the need to explicitly open/read/write files. For example, when writing the forward model algorithm, we don’t really care where e.g. the ion temperatures or gravity vectors come from, be that from a memory cache, disk file, or generated on demand. Three classes, working together with the low-level I/O mechanisms described above, take care of this layer. They are the DataLoader template class, the Generator concept, and the Cache template class. To achieve this, each logical dataset (which could be a vector, a matrix or anything else) is assigned a unique identifier, which, in a user friendly manner, describes the data, and where it may be stored. An example of such an id 106 is gaim physics 20010921 023600.mat%ion temperatures vox, which identifies the ion temperatures vector in voxel centers as of 2001/09/21 2:36 UT, also indicating that it might be saved in the file named gaim physics 20010921 023600.mat. In GAIM++ the data identifier is encapsulated in the class DataId. Note that the mechanism of the Loader is independent of the actual implementation of the identifier (and is templated on it), so a copletely different implementation could be used instead. Given a data identifier, the application code can request to Load it using a DataLoader. The DataLoader has an internal cache (using theCache class), so it first checks if the cache contains the requested ID. If it does, the associated data is returned to the caller. Otherwise, it checks if the data is present in the disk file. Again, if it’s there, it loads it, puts it in the memory cache and returns it to the caller. In case if it’s not in the file either, the DataLoader has an internalGenerator, which as the name indicates is in charge of generating the data in question. Thus, as a last resort, theDataLoader asks the Generator to create the requested data on the fly (one will recognize the Delegation and Builder design patterns here). Once the generator completes that task and returns the data to the loader, it is put in the memory cache, and, depending on the policy set on the generator, it is saved in the disk file for a use at a later time. It’s easy to see how this design can also be laveraged for running and generating some data off-line, then using it later for real-time processing. For example, many of the datasets in the physics input (see section 8.3) can be pre-generated and used later in the forward model. The same way, once the forward model is ran, its results can be used in different data assimilation modules, without the need of generating them internally. 107 This design of the DataLoader, has one more advantage: it allows more com- plex abstractions built on top of it, shielding the application layer even more from the necessity of dealing with messy I/O realted bookkeeping. GAIM++ includes a couple of classes in this category. As an example, let’s examine the InterpolatingDataLoader class. As the name suggests, this loader is used when there is need to interpolate between 2 data sets. For example in the physics mod- ule, the data is generated and saved at 12 minutes intervals only. However, all the other modules will need these physics quantities at arbitrary times. not necessarily aligned with the 12 minutes time grid. This is where the InterpolatingDataLoader comes into play: it hides all the details about implementing the interpolation, and makes it completely transparent to the caller. The caller does not even have to be aware that there is interpolation going on in the back stage. Thus if in the future, our physics module is capable of generating the same data with no interpo- lation needed, the model code will not have to be changed at all (only the loader will need to be changed inside the physics module). The actual implementation of the InterpolatingDataLoader and other such classes follows the Composite design pattern, i.e. it implements the DataLoader interface, and has an internal member of typeDataLoader on its own. Two other classes in this category of data loaders in GAIM++ areLinearDataLoader andPerturbedDataLoader. The added benefits from these layers of abstraction come with a relatively small cost: because of the number of indirections the data goes through, understanding the actual flow of the data or debugging may become more complex, especially when the higher level loaders are nested within each other. To minimize the possibility of inadvertently modifying some cached data, that another part of the application does not expect to be modified, theLoad() method of theDataLoader returns the data 108 only as an immutable object (a design decision, made after a couple of sleepless nights of debugging). 8.1.4 Other utilities GAIM++ comes with many utility class. Among them: Namelist - allows parsing and accessing fortran style namelist files. Cache - provides policy based resource caching. Different policies are imple- mented for cache expiration (such as LFU - least frequently used, LRU - least recently used), other policies are implemented for determining when the cache limit is reached (such as based on memory occupation or the count). Different logging mechanisms can also be plugged into this class. DbgMessage - provides debug, log and error message capabilities, using log/error levels, categorization, etc. Debug statements are 0-cost when GAIM DEBUG is not defined. GaimExceptions - GAIM uses a hierarchy of typed exceptions, derived from std::exception. No other types of exceptions are thrown in the code. Time - provides time related functions, with conversions from/to different for- mats. StringTokenizer - an efficient (but a bit unsafe) string tokenizer imple- mentation, with a pluggable token finder. DelimLayout - provides an abstraction for a delimitted text record format, with methods to access and convert each data field (based on the field name or index). 109 BSpline - a program that can perform 4 different tasks, related to spline func- tions. – Evaluate the values of spline basis functions at some bin points. This is useful for visualizing the spline basis functions – Given coefficients for the basis functions, evaluate the values of the spline function formed with these coefficients at some bin points. This is useful for visualizing spline functions. – Given some points on the plain, calculate the coefficients of a spline func- tion, which would interpolate these points. – Given some points on the plain, calculate the coefficients of a spline func- tion which has the best fit for these points. The program can handle splines of any order, as well as “periodic” splines. Many other utility programs are included in GAIM++ for generating density pro- files, transforming different geometry grids, filtering and simulating observation data, etc. 8.2 Geometry Module As discussed in Section 4.1, GAIM++ uses the finite volume method to represent the partial differential equations which model the Earth’s ionosphere. The space around the Earth encompassing the region of the ionosphere that is modeled, is partitioned into a meshed geometry. The 3-dimensional “finite volumes” thus obtained are called voxels. The geometry module contains classes that model these abstractions. It is 110 in charge of providing all the quantities needed by the physics, forward or assimila- tion modules. Examples of such quantities are the coordinates centers of the voxels and faces, volumes of voxels, areas of faces, distances between different locations, neighbor voxel or face information, etc. We briefly discuss now the design of some components in the geometry module, and the decisions that lead to that design. Coordinate System - GAIM++ uses five different coordinate systems throughout all the modules. They are PQL - representing the Earth’s geomagnetic and geopotential field lines. Spherical or RLL - the spherical coordinate system with Earth’s rotational axis as its vertical axis. Magnetic Spherical or MRLL - the spherical coordinate system with a vertical axis being parallel to Earth’s magnetic axis. Cartesian or XYZ - the standard Cartesian coordinate system, with Earth’s rota- tion axis as the Z axis. Magnetic Cartesian or MXYZ - the standard Cartesian coordinate system, with Z axis parallel to the Earth’s magnetic axis. Note that for PQL, MRLL and MXYZ, depending on the Dipole used (see section 2.3) the origin of the coordinate system may be shifted away from the center of the Earth. TheCoordinateSystem class represents one of these five different coordinate systems. This class is a simple marker class (or traits class), so it does not actually carry real data, or expose an interface. Early in the design of the geometry module, some effort was spent in the direction of implementing a unit/dimension library to be used as part of GAIM++ geometry and physics modules. This task proved to be not an easy 111 one, and the project started stalling because of it. So it has been abandoned. However, some ideas from that effort have found use in the implementation as exemplified by the coordinate system and direction traits. There are several ongoing independent efforts to implement an efficient (0-cost) units/dimensions library (e.g. in boost [8]), and when they reach maturity, they could provide important benefits in developing a geometry or physics modules for numerical models. Point - This class represents a point in 3-dimensional space in one of the given coordinate systems. It provides methods to convert to and from any other coordi- nate system. The first attempt at this class used to include all the coordinates for all 5 systems. This proved to be unpractical and inflexible from several points of view (space efficiency, updating changes in all the systems when one of them is changed, etc). So we decided to split the class into 5 different classes. It was only natural to make the class a template so it can handle any coordinate system (and so we wouldn’t have to write/maintain multiple copies of essentially the same code). The actual code for doing coordinate transformation is delegated to a separate class - CoordinateConvertor, which can be configured to use different types of dipoles (centered, tilted or eccentric). We also decided to allow implicit conversion between the differentPoint classes. Eventhough C++ implicit conversions are generally con- sidered error-prone, in our case it will not lead to errors, since the conversion does not change the logical meaning of the object, it still represents the same logical point in space. The advantage for this approach is that, templetized routines can be written which act on any type of point, and they can still seamlessly call any routines that requires a point in specific coordinate system without unnecessary code clutter. The 112 only drawbacks for implicit conversions in our case could be 1) performance degrada- tion, if there are lots of conversions in some tight loops that are happening unintention- ally. This can easily be fixed by simple profiling after (and only after) the program is proved to work properly, and proved to spending too much time in the conversion rou- tines. 2) numerical errors, since they accumulate with every conversion. Again, once the program is running correctly, it should be examined to make sure that there are no unnecessary conversion back and forth between different coordinate systems. Even though the point class has been split, and contains the representation in one coordinate system only, this design allows a pretty straightforward way to create new classes, which would contain the representation in several systems. The ”Composition” or ”Multiple Inheritance” design patterns can be used to implement such a class. In fact, GAIM++ contains an implementation of a generic class that represents a point in 2 coordinate systems (MCSPoint) which uses multiple inheritance. Voxel - As mentioned above, this class represents the basic unit of volume used in our modules. It contains information about the size of the voxel, its volume, the center point, the neighboring faces and voxels. It also has a flag specifying if the voxel is ”active” in the current geometry (see section 4.1) for discussion about active voxels). Face - This class represents one of the 6 faces of a voxel. It contains information about the size of the face, its area, the center point, the neighboring voxels and the normal vector. TheFace class is templetized on the dimension trait (e.g. P, Q, L). The primitive objects above (Point, Face, Voxel) are generated based on the meshed geometry, as noted above. The geometry module has some classes to represent this mesh (or grid), as well as the actual region of interest that is being modeled. IntervalPartition - This class represents a partition of an interval into seg- ments. It has methods to access the endpoints of the interval and the segments and their 113 centers. It also has methods to efficiently find the interval where a given point is. The intersting thing to note about this class is that an interval partition can be ”periodic”, i.e. this class can also represent the partion of a circle. This is used to represent the lon- gitude direction in our geometry. Like thePoint class above this class is templetized on the dimension trait. Grid - This class represents a rectangular grid in 3-dimensional space. It’s com- posed of 3 interval partitions, one for each dimension. The class is templated on the coordinate system, so it can represent any kind of grid (PQL, RLL, etc). Region - This class represents the actual region in space that is being modeled. In essence it is a collection of V oxels (and their respective faces, etc), that represent the region of interest. The region is initialized based on a grid and some geometry constraints about the part of the ionosphere to be modeled. Given the information, Voxel objects are generated to represent each ”volume element” of the grid. Thus logically, Region is a ”subdomain” of the grid (the grid is rectangular, while the region doesn’t have to be). The class has methods to save or load the entire information about the geometry in specialized binary files (e.g. Matlab orNetCDF files). It also takes care of the proper ordering of the voxels that is suited for the numerical problem at hand. The geometry module has a few other classes, that instead of encapsulating data about geometric objects, encapsulate some mathematical concepts dealing with geom- etry, such as distance, smoothing, interpolation, projection, etc. An important point to note, is that the classes in the geometry module encapsu- late only geometry related information. No physics related data (such as density, or gravity) is encapsulated in this objects. This might seam to be against the concepts of 114 encapsulation from a purely object-oriented point of view. However, from a numeri- cal .... point of view, it is important to keep the data that will be used later in matrix operations in a ”native” format (e.g. vectors, matrices). 8.3 Physics Input Module The Physics module is in charge of providing input data to the Model. This data includes information such as gravity, temperatures, chemical reactions, etc. It is impor- tant to note, that each one of those quantities usually comes from a model on its own. For example the information about the horizontal wind comes from the HWM93 [15] model. Some of the inputs can have alternative models that can be used, such as MSIS-90 [16] or NRLMSISE-00 [26] for the densities and temperatures of the neutral species. Also, clearly as science and experiments advance, each of these areas will come up with better models. So one of the design goals of the physics module was to make it possible to easily swap different model-implementations of each of the phys- ical quantities. Thus, in GAIM++, each of the physical quantities, is represented by a pair of Interface-Implementation classes (e.g. Wind-WindHWM). In addition, these interfaces are wrapped intoDataLoader’s described in section 8.1.3. ThePhysics class itself only works with these loaders, it’s completely unaware of the implementa- tions. So, when an alternative implementation is needed, no code will need to change in the Physics class, or the other Equation and Data Assimilation modules. The use of DataLoader’s also allows the physics model to easily fit into the 4DV AR process, where the physics drivers need to be changed many times during the mini- mization process. This is done by swapping the default (“climate”) DataLoader’s for those drivers with other loader objects (of classPerturbedDataLoader). 115 Figure 8.2 shows more details about each component of the physics module and how they are fed into the forward model. The physics module also includes a standalone program PreparePhysics, to pregenerate physics input data. 8.4 Forward Model The Forward Model module is in charge of “evolving” the ion density state from a given time point to the next one. This is done by constructing all the components in the matrix equation described in section 4.1, then solving the matrix equation thus produced. The implementing class (MultiIonEqution) contains several methods each in charge of filling in some parts of the matrices or vectors involved in (4.12), such as convection velocities, solar production, etc. Figure 8.2 shows the details on the flow of the program. The actual solving of the matrix equation is handled through uBLAS-Lapack binding library (also part of uBLAS), which allows the use of highly optimized/platform specific Lapack implementations for efficiently solving the equa- tions. Typically the forward model is called to evolve the state in fixed size time steps (e.g. 12 minutes). However in general the model cannot be solved for such a time interval in one step (see section 4.1). Therefore the Forward Model internally determines the maximum time step than can be taken at once, then iterates as many times as needed, to arrive to the desired time point. 116 Physics Input Solar Production (V) ● Solar EUV irradiance ➔ EUV94 ● Optical Depth ➔ empirical ●Gravity (V) ➔ empirical ●Neutral Temperatures (V) ●Neutral Densities (V) ➔ MSIS 90 ➔ NRL MSIS 2000 E×B Drift (p) ●E×B Drift (eq) ➔ empirical ➔ FejerScherliess Parallel Velocity (q) (convection) ●Electron Temperature (V,q) ●Ion Temperatures (V,q) ➔ empirical ●IonNeutral Collisions (q) ●IonIon Collisions (q) ➔ empirical ●Gravity (q) ➔ empirical ●Horizontal Wind (q) ➔ HWM93 Chemical Reactions (V) ●Neutral Densities (V) ●Neutral Temperature (V) ➔ MSIS 90 ➔ NRL MSIS 2000 ●Ion Temperatures (V) ➔ empirical Forward Model Flux on Psurfaces Flux on Qsurfaces (explicit component) Production due to Chemical Reactions Production due to Solar Radiation Flux on Qsurfaces (implicit component) Parallel Velocity (q) (diffusion) Geometry Figure 8.2: The flow of the GAIM++ forward model. The subscripts (V),(p),(q) indi- cate whether the corresponding quantity is needed on the center of the voxel, q-face or the p-face, respectively 117 This module is also capable of saving the matrices computed for the state evolution equation into disk files, so they can later be used in the data assimilation modules (e.g. Kalman Filter). 8.5 Observation The data assimilation modules in GAIM++ are required to handle several different types of observation data, such as absolute and relative TEC, in situ densities, etc. These observations may be stored in disk files, or may be streaming. In either case the input can have many different layouts, formats and data content. The Observa- tion module is in charge of handling all these differences, and abstracting out these details from the rest of the application, by providing a simple interface to them. It also provides hooks to make the handling of new types of observations easy. Fig- ure 8.3 shows the Observation class hierarchy. The Observation class is a virtual class: its subclasses are required to implement a few methods, such asRead() or GetObsMatrix(). Based on this, it provides some utility methods to get other observation related data structures (e.g. residuals, innovation, etc). On the next level down, the FileBasedOservation class handles observation data stored in disk files. It is templated on the type of measurements taken. It implements theRead() method above, based on the abstractReadObservationLine() method. Its main task is to take care of details like files sorted backwards in time, or the bookkeeping of an extra record read when determining where to stop the current read operation. But, because it needs theReadObservationLine() method, this class is also abstract. However, the classes which inherit from it, would now only have to implement a simple method to read just one line out of the observation file into the measurement data struc- ture. An example of such a subclass is GpsTecObservation. This is a concrete 118 class, and is specially made to handle the particular format and content of the GPS TEC DUMP files. It uses theDelimLayout class described in 8.1.4 to parse and store the observation data into the GpsTecMeasurement structures. It is also in charge of interpreting that data (since only it has the knowledge about its meaning) and creat- ing the observation matrix based on it (thus implementing the GetObsMatrix() method). The SerializedObservation class shown in figure 8.3 is another example of the Composite pattern: it implements theObservation interface, and internally has an instance of FileBasedObservation. This class is used for saving the observation data needed by the data assimilation programs in a binary format, which makes the repeatedRead() operations very efficient. GpsTecObservation InsituObservation Observation SerializedObservation EuvObservation RelTecObservation FileBasedObservation<M> Figure 8.3: Observation class hierarchy 119 8.6 4DV AR As was described in section 6, the 4DV AR is an iterative process, where the estimated parameter is refined on each iteration. Figure 8.4 shows the logical flow of this iterative algorithm. Solve Adjoint Equation (backwards) Going backwards, for each time step ● read the observations for the time step ● calculate residual and innovation vector g k ● accumulate the norm of the residual in cost function J ● solve adjoint equation using A k -1 and B k to find the co-state v k ● calculate F k j for each element of the parameter α j ● accumulate the sum v k F k j in ∇ j J Δt 1 Δt 2 Δt N Perturb the physics drivers with α Check the cost function. Exit if small enough. Start with initial parameter α=α 0 Solve Forward Equation ● save state n k and matrices A k and B k on each step ● mark the time steps Δt 1 Δt 2 Δt N Obs. Data Obs. Data Obs. Data Run L-BFGS iteration to update α Figure 8.4: The flow of one cycle of the 4DV AR process 120 Even though the ionosphere has the natural period of one day and our forward solver is typically ran for 24 hours, because of the very dynamic nature of the iono- sphere and the underlying physics forces, only one value of the parameter vector for the entire day cannot be used to adjust the drivers to a satisfactory degree (in fact the process may not even converge if ran for long periods of time at once). Thus, GAIM++ uses shorter time slices, called cycles, to run the iteration. For each cycle the process uses the observations falling into that time period and finds the minimizing parameter. It is important to note, that the length of a cycle has a lower boundary as well, dic- tated by the amount of observations available. This is because if there are no enough observations for the given cycle, the process may not converge. Different strategies for slicing the day into cycles can be used in GAIM++, including overlapping cycles. The rational for using overlapping cycles, is that the reuse of the observation data falling into the overlap region of time, makes the transition between the solutions from one cycle to the next one much smoother. Typical cycle configurations used in GAIM++ are: 2 hour cycles with no overlap, 3 hour cycles with no overlap, 4 hour cycles, with 1 hour overlap. Another important aspect of using the cycles, is that each cycle uses the minimizing parameter of the previous cycle as its initial value. This helps reducing the number of iterations needed to converge, since the solution of the previous cycle is likely to be close to the solution of the next cycle. As can be seen from figure 8.4, after each iteration, the algorithm checks the con- vergence criteria, which consists of several checks. Namely, the iteration is stopped if i) the cost function value itself is smaller than a predefined threshold (adjusted by the number of observations used for the cycle); or 121 ii) the decrease in cost function is smaller than a predefined small fraction of the first cost function value. This is because it is considered that the first cost func- tion value should be a “reasonable” value already, since its corresponding state is the result of a physics based model; or iii) the gradient of the cost function is smaller than a predefined small fraction of the cost function value (this check is only done in certain situations related to the specifics of the L-BFGS algorithm), or iv) the maximum number of allowed iterations is reached. In figure 8.4 the two central blocks are responsible for one evaluation of the cost function. When solving forward, the states and the corresponding matrices are saved using the mechanisms described in sections 8.1.2 and 8.1.3. This provides both effi- ciency and accuracy in solving the adjoint equation compared to the alternative of recomputing these matrices on the fly. When going backwards, the Observation class from section 8.5 is used to read in the observation data, and compute the residual and innovation vectors. As mentioned before in 8.5, since the same observation data is read for each iteration (within a cycle), the original observation files are serialized into an efficient binary format and saved at the very beginning of the program, before the first cycle starts. The adjoint equation is solved using Lapack’s banded solver rou- tine gbtrs (connected to uBLAS through its lapack bindings). Then after, the cost function gradient is calculated (more about this below), the results are fed into the L-BFGS routine (rewritten in C++ by scitbx [24]), which is in charge of supplying the next candidate for the parameter. After the convergence check described above is done, the algorithm loop is restarted, if needed. Before the forward solver can pro- ceed again, the drivers have to be adjusted with the new parameter. The central class 122 which coordinates the handling of the parameters, adjusting the drivers and calculating the components needed for computing the cost function gradient isPerturbation. Figure 8.5 shows the Perturbation class hierarchy and its relationship with the PerturbedDataLoader class. Perturbation is a virtual base class, it leaves the job of implementing the main functions (such asGetPartialA(), in charge of calculating @A @ j from section 6) to its subclasses. However, based on these functions, it does implement theGetF(), function in charge of computingF jk (see section 6). And sinceF jk is the only thing needed by the cost function evaluation function, that function is made to interface with the abstractPerturbation class only, effectively shielding it out from the details on how the actual computation of @A @ j , @B @ j , @P @ j and F jk are done. ProductionPerturbation WindPerturbation Perturbation PerturbedDataLoader DriftPerturbation MultiPerturbation PhysicsPerturbation AbstractPhysicsPerturbation DataLoader Gaussian3DPerturbation Physics Figure 8.5: ThePerturbation class hierarchy 123 Another benefit of this abstraction is that using theMultiPerturbation class, without any change to the main 4DV AR loop, we are able to handle adjustment to multiple drivers. Internally, this class does all the bookkeeping (like index remapping) needed for handling multiple perturbations, but publicly, it provides the same sim- ple interface ofPerturbation. Note thatMultiPerturbation is yet another example of the Composite pattern: it contains a list of other Perturbations while implementing thePerturbation interface itself (looking even more closely, MultiPerturbation is also an example of CRTP (Curiously Recurring Template Pattern), but this is more of a side effect, than a design need). Back to the 4DV AR cycle, when the minimization routine supplies the new candi- date parameter, it is set into thePerturbation object. As described above in 8.3, the Physics class uses DataLoader objects for supplying the data related to the drivers we are adjusting. After getting the new parameter,Perturbation swaps out these loaders and replaces them by its own instances ofPreturbedDataLoader. Thus, later, when the call comes down from the forward solver through the modified Physics object to thesePreturbedDataLoader’s, they will load (or rather gen- erate, see 8.1.3), the new physics quantities based on the new parameter. Note that, PreturbedDataLoader is a wrapper class, which delegates the actual job of per- turbing the data back to thePerturbation class (it can be thought of as an example of the Proxy pattern). The main remaining difficulty in this loop (and, actually, in the 4DV AR process in general) is the evaluation of the partial derivatives @A @ j , @B @ j and @P @ j (see sec- tion 6). Since the computation of each of these quantities needs intimate details about the corresponding driver and its specific perturbation, this job is that in the concrete classes ProductionPerturbation, DriftPerturbation and 124 WindPerturbation. Moreover, these classes need some extra information from the forward solver, which it (very helpfully) pre-computes and saves during the for- ward run. Then, these classes use the formulae described in section 6.2.1 for their implementions. Note that, since in our implementation both the Production and Wind perturbation use a Gaussian 3D perturbation (see section 6.2.1), they inherit from a base class, which has a generic implementation of these types of perturbations. The classGaussian3DPerturbation makes use of theRegionTransformer class from the geometry module to do the Gaussian interpolation (6.13). 8.7 Kalman Filter The Kalman Filter ([20]) is an alternative method of data assimilation implemented in GAIM++. Instead of adjusting the physics drivers as 4DV AR does, it directly adjusts the density state based on the observation data, by balancing it against the state obtained from the forward model. This balancing is controlled by the observation noise, process noise and the state covariance matrices. One of the design goals of the Kalman Filter module in GAIM++ was to sup- port different types of state covariance matrices as a way to control the computational intensity of the process. Among them diagonal, banded, sparse and full. GAIM++ does this by making the main class for the Kalman filter a template class, templated on the type of the covariance matrix. Since the uBLAS library has support for all the different matrix types that we need to handle, and they all present more or less the same interface (as they are all models of the matrix concept), it makes this task easily manageable. It is important to note however, that because the Kalman Filter loop involves an update to the covariance matrix both in the time propagation of the covariance and in the Kalman update step (see figure 8.6), it is strictly speaking not 125 possible to maintain a non-full covariance matrix. So in practice, GAIM++ imple- ments a so-called “structure-preserving” Kalman Filter (referring to the structure of the covariance matrix). It basically consists of throwing away all the elements of the updated covariance matrix that were not part of its structure before the update. Thus a diagonal covariance matrix will always stay diagonal, a banded one will stay banded with the exact same diagonals, etc. In reality of course, instead of computing all these throw-away elements then discarding them, GAIM++ simply skips their computation. For that purpose, GAIM++ has a set of matrix manipulation routines (mainly versions of the well known axpy routine) that are specialized for the structure-preservation. Also note, that the structure-preservation is not necessary (and is not done) when run- ning with a full covariance. Thus the non-optimality of the Kalman Filter introduced by these technique is just a trade off between computational efficiency and the correct- less of the filter. Another specific aspect of our implementation is the handling of the three covari- ance matrices involved in the Kalman filter process, that is the state covariance matrix we discussed above, the process noise covariance matrix and the observation covari- ance matrix. State covariance matrix It is often said that the art of implementing a successful Kalman filter is in provid- ing a good state covariance matrix. GAIM++ allows the user to supply an initial covariance matrix to the process. This is useful for testing purposes, for restarting a previously interrupted process, as well as for using a covariance matrix obtained by external means. When none provided however, the program internally constructs this initial matrix, based on the pre-specified correlation lengths. The correlation lengths 126 are distance measures which specify the “strength” of the relationship between vox- els as a function of their distance from each other. Typically this strength decreases exponentially with the distance (or, more precisely with the square of the distance). So, using the correlation lengths the program first constructs the correlation matrix. Then the covariance matrix is constructing by making each value proportional to the corresponding correlation and to the ion densities in both voxels involved. The intu- ition behind this choice, is that the uncertainty is higher for voxels with high density. To allow some meaningful uncertainty for voxels with very small densities, a small additive term is added to the diagonal of the covariance matrix. Because of the big variability in voxel sizes in our grid (see section 4.1) we introduce the so-called voxel weighting in our covariance calculation. The motivation behind this is as follows. The ion density we assign to each voxel represents the average of the real ion density in the entire voxel. Thus, the bigger the voxel, the less uncertain this average should be. Therefore we make the covariance to also be inverse proportional to the average of the volumes of the two voxels involved. Process noise, aka “Q-bump” The process noise term used in GAIM++ follows the same algorithm described in the above discussion about the initial state covariance matrix. In fact in our current implementation these two are one and the same, with the initial covariance being equal to theQ-bump using the initial ion density state. Observation noise The observation noise includes the measurement errors and the representation error. GAIM++ Observation class provides an interface for both of these errors. This 127 allows the implementing subclasses to customize the algorithms so that they are suit- able for the particular type of observation being implementing. For example, currently in our implementation of GpsTecObservation described above, we use a pre- configured constants values for both the measurement error and representation noise. An alternative implementation for the representation noise is being investigated how- ever, where the representation noise depends on the observation matrix and the current density state. Namely, letR be the representation noise (diagonal) matrix, letH be the observation matrix andn the ion density state. Then we set R(i;i) = X k (H(i;k)(n(k) +)) 2 wherei andk represent the observation and the state vector indices respectively, and and are some constant. Without going into much details, the motivation behind this formula, is that the representation noise should depend on the way the GPS observation link intersects our discrete grid. Therefore it should be a function of the density portion contributed to the link from that voxel (represented in the formula above byn(k)) and the so-called partials - the segments cut from the link by our 3D voxels (represented in the formula above byH(i;k)). All of the noise terms discussed above are subjects of active research and bound to be changed and finetuned as a result of it. 128 Read observation batches N max N max N max N max N max N max Going forward, for each time step ● read a batch of maximum N max observations in the time step ● calculate the observation matrix H and residual dy ● get the measurement noise covariance R ● calculate the Kalman update matrix: K = PH t (HPH t +R) -1 ● calculate the Kalman update to the state and the covariance n = n + Kdy and P = P - KHP ● constraint the state by replacing negatives with 0. Solve Forward Equation ● save state n k and matrices A k and B k on each step ● record the time steps the forward solver took ● Propagate covariance P in time using A k and B k P = (A k -1 B k ) P (A k -1 B k ) t ● Add process noise to the covariance: P = P + Q k Δt 1 Δt 2 Δt N Δt 1 Δt 2 Δt N { { { { { { Figure 8.6: The flow of the Kalman filter process Figure 8.6 shows the logical flow of the Kalman Filter algorithm as implemented in GAIM++. As in the case of 4DV AR, the forward solver saves the state transition matrices when evolving the density state, so that they can later be used for the covari- ance matrix propagation. Also similar to the 4DV AR module, the Observation 129 class is in charge of reading in the observation data, and supplying the residual vec- tor needed in the algorithm. Because the dimensions of the matrices involved in the Kalman update step are proportional to the number of observations read, to have con- trol over the computational intensity of that step we need to be able to control how many observation records are read before the Kalman update is invoked. Helpfully, the Read() function of the Observation class can limit the number of observa- tions read at a time. Even though this technique solves the practical problem of solving large linear systems, it introduces another slight non-optimality in the Kalman Filter process itself, in the case where a non-full covariance matrix is used. As described above, this is because after each batch of observations read, a part of the covariance matrix is effectively thrown away in the Kalman update step by the process described above. When propagating the covariance matrixP in time (see figure 8.6), the inverse of A needs to be computed. This matrix could be rather large, since its dimensions are the size of our density state. Therefore the computation of its inverse could be expensive in both time and memory. We solve this issue by recalling that A is a diagonally dominant (tri-diagonal) matrix with the diagonal being significantly larger than the rest of the elements. For this kinds of matrices,A 1 can be approximating with a good 130 accuracy byD 1 D 1 (AD)D 1 , whereD is the diagonal ofA. This follows from the following approximation, by noting thatD 1 (AD) is small A 1 = (D + (AD)) 1 = D 1 +D 1 (AD) 1 = 1 +D 1 (AD) 1 D 1 = 1D 1 (AD) + D 1 (AD) 2 D 1 1D 1 (AD) D 1 = D 1 D 1 (AD)D 1 With this approximation, the computation of A 1 is very efficient, and results in a tri-diagonal matrix, which later makes computingA 1 B efficient as well. 131 Chapter 9 Conclusion As a continuation of the efforts in developing an accurate global ionospheric mon- itoring and forecast model, the work presented in this thesis represent a significant enhancement of the earlier version of GAIM. By addition all most relevant ion species in this version of forward model we can examine the interactive processes occurring in the plasma of the Earth’s ionosphere. Our results have shown that the numerical approach taken in approximating the solutions of the ionospheric model equations leads to physically meaningful and stable ion densities. The model compares favorably to other existing numerical models of the ionosphere. Moreover, like our earlier version of the forward model, this new version of the model has the usual features that are particularly useful for data assimilation. In particular, the parameterization of the driving forces and the solution to the adjoint equation allow efficient implementation of the 4DV AR approach, as well as, Kalman filter approach for data assimilation. Our work on 4DV AR demonstrates that 4DV AR approach can be effectively used to estimate ionospheric driving forces. However, we have also shown that when our observation is limited to ground based slant TEC, there are limits on what driving forces can be unambiguously resolved. This is entirely to be expected for a complex dynamical system like the Earth’s ionosphere. The development of new ionospheric model and 4DV AR tools also lead to opportu- nities to further investigation approaches for more accurate ionospheric monitoring and 132 forecast. There are many outstanding issues with the basic modeling of ionosphere, as well as, the assimilation of data. From a modeling point of view, the following areas need to be further explored. Change of models for the Earth’s basic magnetic field. The inaccuracy of a dipole model causes systematic errors in the location of some of important fea- tures of the ionosphere including equatorial anomaly. The direct use of IGRF can significantly reduces this error. Creation of geometry more closely resembles spherical grid. The horizontal scale size of the Earth’s ionosphere is substantially larger than the vertical scale size. A more spherical like modeling grid allows better balance between model resolution and the overall dimension of the model. Combination of regional and global models to allow high spatial resolution model over region of high interest. Improvement in the numerical approaches in solving the coupled nonlinear con- vection and diffusion equations to allow large integration step sizes. On the 4DV AR side, our work represents the beginning of investigation of the potential of this approach for data assimilation. Unlike in an OSSE, in an actual use of 4DV AR with real data, there is no independent measurements of the driving forces. An alternative metric of merit must be developed to determine whether or not the estimated driving forces are more accurate than climatological values. In fact, the greatest promise of 4DV AR is to improve the model’s prediction skill by estimating the drivers. Considerable experience with 4DV AR must be gained in tuning and calibrating the model to produce reliable forecast. 133 References [1] D. Anderson, A. Anghel, K. Yumoto, M. Ishitsuka, and E. Kudeki, Estimat- ing daytime vertical exb drift velocities in the equatorial f-region using ground- based magnetometer observations, GEOPHYSICAL RESEARCH LETTERS 29 (2002), no. 12. [2] D.N. Anderson, A theoretical study of the ionoshperic f region equatorial anomaly–i. theory, Planetary Space Science 21 (1973), 409–419. [3] , A theoretical study of the ionoshperic f region equatorial anomaly–ii. results in the american and asian sectors, Planetary Space Science 21 (1973), 421–442. [4] G.J. Bailey and R. Sellek, A mathematical model of the earth’s plasmasphere and its application in a study of he(+) at l = 3, Annales Geophysicae 8 (1990), 171–190. [5] G.J. Bailey, R. Sellek, and Y . Rippeth, A modelling study of the equatorial topside ionosphere, Annales Geophysicae 11 (1993), 263–272. [6] P.M. Banks and G. Kockarts, Aeronomy, part b, Academic Press, New York and London, 1973. [7] S. Benharbit, A. Chalabi, and J.P. Vila, Numerical viscosity and convergence of finite volume methods for conservation laws with boundary conditions, SIAM J. Numer. Anal. 32 (1995), no. 3, 775–796. [8] Boost, Boost c++ libraries,http://www.boost.org. [9] C. G. Broyden, R. Fletcher, D. Goldfarb, and D. F. Shanno, L-bfgs optimization method,http://en.wikipedia.org/wiki/L-BFGS, 1970. [10] J.M. Burgers, Flow equations for composite gases, Applied Mathematics and Mechanics 11 (1969). [11] Z. Cai and S. McCormick, On the accuarcy of the finite volume element method for diffusion equations on composite grids, SIAM J. Numer. Anal. 27 (1990), no. 3, 636–655. [12] G. Colvin and B. Dawes, shared ptr c++ library, http://www.boost. org/libs/smart_ptr/shared_ptr.htm, 1998. 134 [13] R. Daley, Atmospheric data analysis, Cambridge University Press, 1991. [14] G. Hajj, B. Wilson, C. Wang, X. Pi, and G. I. Rosen, Data assimilation of ground gps total electron content into a physics-based ionospheric model by use of the kalman filter, Radio Science 39 (2004). [15] A. E. Hedin, Horizontal wind model, http://nssdc.gsfc.nasa.gov/ space/model/atmos/hwm.html, 1993. [16] A.E. Hedin, Msise-90 model 1990,http://modelweb.gsfc.nasa.gov/ atmos/msise.html, 1990. [17] J.D. Huba, Sami2,http://wwwppd.nrl.navy.mil/sami2-OSP, 2001. [18] F. John, Partial differential equations, Springer-Verlag, New York Heidelberg Berlin, 1982. [19] M.C. Kelley, The earth’s ionosphere: Plasma physics and electrodynamics, Inter- national Geophysics Series 43 (1989). [20] D. G. Luenberger, Optimization by vector space methods, Wiley-Interscience, 1969. [21] L. Mandrake, B. Wilson, C. Wang, G. Hajj, A. Mannucci, and X. Pi, A per- formance evaluation of the operational jet propulsion laboratory/university of southern california global assimilation ionospheric model (jpl/usc gaim), J. Geo- phys. Res. 110 (2005), 10. [22] MathWorks, Matlab,http://www.mathworks.com/. [23] N. Olsen, Terence J. Sabaka, and Lars Toffner-Clausen, Determination of the igrf 2000 model, Earth Planets Space 52 (2000), 1175–1182. [24] PHENIX, Computational crystallography toolbox, http://cctbx. sourceforge.net/, 2002. [25] X. Pi, C. Wang, G. A. Hajj, G. Rosen, B. D. Wilson, and G. J. Bailey, Estimation of eb drift using a global assimilative ionospheric model: An observation system simulation experiment, Journal Geophysical Research 108 (2003). [26] M. Picone, A.E. Hedin, and D. Drob, Nrlmsise-00 model 2001, http:// modelweb.gsfc.nasa.gov/atmos/nrlmsise00.html, 2001. [27] J. B. Rawlings and J. G. Ekerdt, Gnu octave, http://www.gnu.org/ software/octave/. 135 [28] A. Stepanov, C++ standard template library, http://www.sgi.com/ tech/stl/. [29] O. Talagrand and P. Courtier, Variational assimilation of meteorological obser- vations with the adjoint vorticity equation. I: Theory, Quarterly Journal of the Royal Meteorological Society 113 (1987), 1311–1328. [30] Unidata, Netcdf, http://www.unidata.ucar.edu/software/ netcdf/. [31] J. Walter and M. Koch, ublas basic linear algebra library, http://www. boost.org/libs/numeric/ublas/doc/, 2000. [32] C. Wang, G. Hajj, X. Pi, G. I. Rosen, and B. Wilson, Development of the global assimilative ionospheric model, Radio Science 39 (2004). 136

Abstract (if available)

Abstract The Earth's ionosphere plays an important role in wireless communication. One of the primary goals of this research is to improve the techniques for monitoring and forecasting of the conditions of the Earth's ionosphere. The work in this thesis focuses on modeling of the ion and electron density changes in the Earth's ionosphere and the assimilation of remote sensing measurements into the numerical model. The governing equations of the ionosphere consist of equations of fluid for chemically active plasma in the ionosphere. These equations form a system of coupled nonlinear convection-diffusion equations. In this thesis, we establish the well-posedness of the equations in some special cases. We also establish the numerical stability and consistency of the approximation. The resulting numerical model of the ionosphere produces stable and physically meaningful predictions of densities for all major ion species in the ionosphere. The numerical method is used in this research as a basis for the development of a4-dimensional variational (4DVAR) data assimilation approach to help in the estimation of the driving forces in the ionosphere. The 4DVAR technique uses the adjoint approach in the computation of the gradient of performance functional.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Empirical approach for estimating the ExB velocity from VTEC map

PDF

Covariance modeling and estimation for distributed parameter systems and their approximations

PDF

A fully discrete approach for estimating local volatility in a generalized Black-Scholes setting

PDF

Mathematical properties of ensemble Kalman filter

PDF

Finite dimensional approximation and convergence in the estimation of the distribution of, and input to, random abstract parabolic systems with application to the deconvolution of blood/breath al...

PDF

Obtaining breath alcohol concentration from transdermal alcohol concentration using Bayesian approaches

PDF

Population modeling and Bayesian estimation for the deconvolution of blood alcohol concentration from transdermal alcohol biosensor data

PDF

Model, identification & analysis of complex stochastic systems: applications in stochastic partial differential equations and multiscale mechanics

PDF

Estimation of random input to semi-linear abstract parabolic systems with application to quantitative description of drinking behavior based on trans-dermal alcohol concentration

PDF

A general approach to using problem instance data for model refinement in constraint satisfaction problems

PDF

Validation of an alternative neural decision tree

PDF

Prohorov Metric-Based Nonparametric Estimation of the Distribution of Random Parameters in Abstract Parabolic Systems with Application to the Transdermal Transport of Alcohol

PDF

CI approach to spectra of quantum dots

PDF

A Bayesian approach for estimating breath from transdermal alcohol concentration

PDF

Optimal debt allocation using a dynamic programming approach

PDF

Functional based multi-level flexible models for multivariate longitudinal data

PDF

An abstract hyperbolic population model for the transdermal transport of ethanol in humans: estimating the distribution of random parameters and the deconvolution of breath alcohol concentration

PDF

M-estimation and non-parametric estimation of a random diffusion equation-based population model for the transdermal transport of ethanol: deconvolution and uncertainty quantification

PDF

A comparison of classical methods and second order latent growth models for longitudinal data analysis

PDF

Inference for stochastic models of molecular data

Asset Metadata

Creator Akopian, Vardan (author)

Core Title Modeling of Earth's ionosphere and variational approach for data assimilation

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Applied Mathematics

Publication Date 05/06/2008

Defense Date 12/06/2007

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag 4DVAR,GAIM,ionosphere,Modeling,OAI-PMH Harvest,variational approach for data assimilation

Language English

Advisor Wang, Chunming (committee chair), Newton, Paul K. (committee member), Rosen, Gary (committee member)

Creator Email vakopian@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m1232

Unique identifier UC1285463

Identifier etd-Akopian-20080506 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-71573 (legacy record id),usctheses-m1232 (legacy record id)

Legacy Identifier etd-Akopian-20080506.pdf

Dmrecord 71573

Document Type Dissertation

Rights Akopian, Vardan

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu