Optimizing the Retrieval of Molecular Concentrations from the Remote Measurement of Flare Emission

Paper from the AFRC 2013 conference titled Optimizing the Retrieval of Molecular Concentrations from the Remote Measurement of Flare Emission by R. Panfili.

Petroleum refineries utilize open-flame combustion of hydrocarbons by industrial flares to eliminate waste gases which are not commercially viable. Incomplete combustion will result in the emission of trace amounts of ozone-forming highly-reactive volatile organic compounds (VOCs) and human carcinogens. Minimizing emissions of VOCs and human carcinogens from flare emissions is driving the need for continuous and autonomous monitoring and combustion flare control technology. Monitoring flare emission is challenging due to the low concentrations of the species of interest and varying environmental conditions. Properly interpreted hyperspectral imagery provides a method to achieve these goals without interfering in the operation of combustion flare facilities. We discuss the development of a retrieval algorithm as a means of determining molecular column amounts from observed radiance with the intent of using those column amounts to compute the efficiency of combustion. The algorithm is based on a non-linear least-square iterative fit to the radiation transport equation.

*[e-mail]; [phone]; spectral.com Optimizing the Retrieval of Molecular Concentrations from the Remote Measurement of Flare Emission R. Panfili*, P. Vujkovic-Cvijin, and H. Dothe Spectral Sciences, Inc., 4 Fourth Avenue, Burlington, MA 01803 ABSTRACT Petroleum refineries utilize open-flame combustion of hydrocarbons by industrial flares to eliminate waste gases which are not commercially viable. Incomplete combustion will result in the emission of trace amounts of ozone-forming highly-reactive volatile organic compounds (VOCs) and human carcinogens. Minimizing emissions of VOCs and human carcinogens from flare emissions is driving the need for continuous and autonomous monitoring and combustion flare control technology. Monitoring flare emission is challenging due to the low concentrations of the species of interest and varying environmental conditions. Properly interpreted hyperspectral imagery provides a method to achieve these goals without interfering in the operation of combustion flare facilities. We discuss the development of a retrieval algorithm as a means of determining molecular column amounts from observed radiance with the intent of using those column amounts to compute the efficiency of combustion. The algorithm is based on a non-linear least-square iterative fit to the radiation transport equation. Keywords: Industrial Pollution, Remote Monitoring 1. INTRODUCTION The process of refining oil produces large quantities of natural gas, presenting a serious explosion hazard. This hazard is mitigated by on-site, open-air burning of the natural gas, a process referred to as flaring. For environmental and health reasons, it is important that flare combustion be as efficient as possible. Incomplete combustion of natural gas results in the release of hazardous organic compounds, which pose serious risks to the environment and to public health [Seebold, 1997]. Examples include a number of known human carcinogens (such as formaldehyde) as well as a large number of trace species which significantly degrade air quality by elevating local ozone levels. These trace species are commonly referred to as highly reactive volatile organic compounds (HR-VOCs), including 1,3-butadiene, butane isomers, ethylene, and propylene. The key to minimizing these adverse environmental consequences is to insure that flares continuously operate near optimal combustion efficiency (CE). Current industry practices do not directly address this issue. Instead, operational conditions are continually monitored based on infrequent visual observations aimed primarily at making sure that the flame is not producing soot (black smoke). Meanwhile, direct measurements of CE using active sensing or gas sampling are cost prohibitive, invasive in that they may interfere with plant operation, and are not continually performed. We aim to address the need to develop a real-time, autonomous, and continuous means of measuring flare CE. The United States Department of Energy has funded research into the development of a real-time, passive, autonomous remote sensing system employing infrared emission spectroscopy to quantitatively monitor the carbon containing flare combustion species and report Flare CE is a key parameter we are seeking to determine. It is often defined by the expression: , 2 2 CO CO HC CO CE        (1) where χ is species mole fraction (or equivalently the volume mixing ratio) in the combusted gas flow (not the active flame), and HC refers to un-burned hydrocarbons. A necessary step to accurately determining CE is the identification and quantification of HR-VOCs that contribute to HC; CE is known to be strongly correlated with HR-VOC production [Seebold et al., 2004, 2007]. We aim to utilize remote detection to ascertain the CE and will illustrate the retrieval algorithm we use in this paper. The use of remote detection for this purpose immediately presents some challenges. First, the fundamental measurement quantity is radiance - the intensity of photons reaching a detector as a function of the frequency of light and time interval. Molecular concentrations cannot be directly measured using passive spectroscopic techniques: they must be inferred from the path radiance. Next, radiance along a single line-of-sight will allow for a measurement of path column amount - the relative concentration of individual molecules along the path - and not molecular concentrations. Finally, CE computations will be dominated by the simplest combustion product hydrocarbons which can be created, namely carbon monoxide and carbon dioxide. Both molecules are present in the ambient atmosphere and atmospheric contributions must be accounted for or mitigated in order to arrive at an accurate measure of CE. These challenges are not insurmountable. First, molecular emission (of infrared light) is a well-studied and well-understood field. Our approach uses the non-linear least-square fitting Levenberg-Marquardt algorithm to derive molecular properties from path radiance. Next, column amounts can be used to compute CE as defined in Equation 1, albeit with values partially consisting of contributions from the ambient atmosphere. Finally, we have investigated approaches to removing atmospheric contributions in an attempt to eliminate its presence in our CE. Our discussion is separated into eight sections. Following this introduction, section 2 provides a discussion of the retrieval algorithm we have used. Next, section 3 outlines the components which make up the matrix of first derivatives (or Jacobian) necessary for our retrieval algorithm. Section 4 discusses our approach towards defining an instrument-specific optimal data set on which our retrieval algorithm will be used. Section 5 provides a demonstration of the retrieval algorithm operating on the instrument-specific optimal spectral data set. Finally, section 6 gives conclusions and future directions for this work. 2. RETRIEVAL ALGORITHM The core challenge of remove detection of CE is to measure an observed radiance and, with that information, determine the conditions from which they arose. To do so requires an inversion (or retrieval) algorithm. Our approach has been to use the non-linear least-square fitting algorithm. Retrievals seek to obtain species densities and temperatures along the radiance path, by minimizing the residual of the observed variable (Fvec):       2 , , = , , , 2 2  i CO CO  comp i obs i CO CO F    L   L    (2) Here, L is the predicted path radiance for a set of guessed conditions while Lobs is the observed path radiance. The algorithm will require a previously-developed forward model [Panfili et al., 2011a] to compute a predicted radiance as well as the derivatives of the radiative transfer equation with respect to all dependent variables. Several algorithms exist to perform this minimization with the Gauss-Newton method [Hartley, 1961] being among the oldest. The Levenberg-Marquardt method [Levenberg, 1944; Marquardt, 1963], which we uses principles similar to Gauss-Newton but can be expected to be better behaved for instances in which an initial guess is sub-optimal. Both methods require only first derivatives in their analysis, avoiding the need to determine higher-order derivatives. They also do not rely on a local linearity of the solution, a basic assumption made, for example, by methods based on series expansion. The Levenberg-Marquardt further possesses advantages over Gauss-Newton for the current problem at hand. Specifically, it is more likely to converge to an accurate solution for a poor initial guess than the Gauss-Newton method. This comes at a cost of slightly longer convergence time. Given the context of this problem, specifically that combustion emission products may span a large parameter range, trading a little extra computation time for greater robustness represents a sensible design decision. This, and other retrieval algorithms, requires derivatives of the equations being optimized. The key derivatives which we will need to construct are that of the derivative of the path radiance with respect to the dependent variables (temperatures and molecular concentrations along the line-of-sight). These path radiance derivative values can be analytically computed for line-by-line algorithms, including the model we have constructed to investigate flare emission [Panfili et al., 2011a]. Our algorithm performs a line-by-line analytical differentiation at runtime with the use of a pre-computed table of temperature derivatives of line strengths. The differentiation of this one term is computed at the same time that cross-sections are generated for our forward model. While functions for path radiance and its Jacobian needed to be created, the core numeric of the Levenberg-Marquardt algorithm did not. There exist robust implementations in the public domain which can be levered and we have made use of one of these implementations. The MINPACK collection [More et al., 1980] provides an implementation of the Levenberg-Marquardt algorithm as well as a series of six interfaces, each optimized for a specific scenario. We are operating with a known function (the radiation transport equation) and we can analytically compute radiance fluctuations, meaning we are able to provide MINPACK with a Jacobian. Furthermore, the volume of data we work with at any one time we are not memory constrained: we will make the tradeoff of increased performance at the cost of slightly more memory. Finally, we do not require flexibility in terms of internal parameter sensitivity as defined in the MINPACK documentation. That is to say, the default settings for certain internal parameters are acceptable. This leads us to conclude that, of the six available interfaces, LMDER1 is the appropriate choice for our needs. MINPACK provides the machinery, but the construction of evaluation routines for both the radiance and derivative values are left to the user. This subroutine is listed as FCN in the MINPACK documentation but can take on any name the user deems appropriate as it is passed as an argument through the interface. FCN has two important roles and control is determined by an integer flag which is also passed through the interface. When the flag is set to 1, the subroutine must return the difference between the computed and observed function values. The results are placed in the fvec array while leaving the Jacobian unchanged.       j obs j comp vec j F  = L   L  (3) When the flag is set to 2, the subroutine must return the Jacobian while leaving the current iterative fit unchanged. This is defined as the derivative of the computed function with respect to the dependent variables.       j comp i j vec i jac i j d dL d F dF      , = = (4) Notice a nuance in these two expressions. The order of differencing the calculated and observed values of the function (in Fvec) is immaterial in computing the least-squares residual. It is, however, significant when computing the Jacobian. We have chosen to subtract the experimental value from the computed value in our implementation. If the order is reversed, the Jacobian would have identical magnitudes, but with a sign flip. The order we chose avoids this unnecessary complication. 3. COMPUTING AND JACOBIAN We will begin by validating the internal radiance computations. Our radiance algorithm assumes the path along our LOS is composed of a series of N homogenous layers. The total path radiance at a single monochromatic frequency can then be expressed as a summation over segment radiance values (R).                  k   path k path k N k k seg k k k seg k N k k N k L N R B t t B t t , 1 , =1 , 1 =1 , =1 =1 1                         (5) In the expression above B is the Planck function and t is the transmission of either a segment along the LOS or the entire LOS, depending on the subscript. Here, we have taken advantage of the following relation between the path transmittance and the segment transmittance, valid for line-by-line methods.   k  seg k k k path t t    , =1 , =   (6) Finally, we note that the transmission through a segment can be expressed as a function of the optical depth (). k   k    k  k   seg t = exp  = exp   (7) Here, the optical depth is expressed as a product of the absorption cross-section () and the column amount (), a product of column density and path length.  =   L =  (8) Three quantities in this expression are computed in the process of finding the path radiance: path transmission, partial path radiance, and segment radiance. They could therefore be stored in memory, awaiting a call to update the Jacobian. This leaves us with the need to differentiate the optical depth and the Black Body function. We begin with the simpler optical-depth term: We cannot analytically solve for the temperature derivative of the cross-section, but we can build an exact value for that derivative at the time our database of cross-sections is created. We will then immediately express the two derivatives we'll need for the Jacobian as: dT d dT d     =  1  (9)     1 =  d d (10) Here, T is the temperature of the segment while each molecular species has its own differentiation with respect to . The key quantities can be expressed in terms of the temperature and column amount derivatives, which we show for the total path radiance below (the Jacobian will be a matrix of similar terms, one for each segment along the line-of-sight).                             m m m m N m m path m m m m m m N dT R B t L L d dT dB dT B dL   = 1        1 (11)                 m m N m m path m m N B t L L d dL    =      1 (12) There remain only two terms which need to be differentiated to complete our Jacobian. The first is the derivative of the black body function with respect to temperature. The second is the derivative of the cross-section with respect to the temperature. The temperature derivative of the blackbody function can be trivially computed from the expression for the Planck function while the blackbody function has no dependence on molecular density.   exp /  1 1 = exp / 2 2 2 2 c T  c T T c dT dB B    (13) Here, c2 is the second radiation constant. The derivative of the cross-section with respect to layer temperature is more involved. The cross-section itself is the sum over all spectral lines (ℓ) and is a product of the cross section for the line-dependent species, the temperature-dependent line center strength (S) and the temperature-dependent line shape profile (f).   ,T  =    S T  f  ,T         (14) The lines used to compute our cross-section originate from either HITRAN [Rothman et al., 2008] or HITEMP [Rothman et al., 2010]. Both databases use the same line notation. They define a procedure to compute line center strength.           2 0 0 2 0 0 2 0 1 exp / = exp 1 1 1 exp / c T c T T T c E Q T S T Q T                           (15) Here, Q is the partition function, T0 is a reference temperature, E↓ is the lower state energy and ν0 is the line center transition frequency. We utilize the Voigt line shape function [Armstrong, 1967], a convolution of the Lorentz line shape which accounts for pressure broadening with the Doppler line shape which accounts for temperature broadening. At the temperatures, pressures and frequencies being considered, the Voigt line shape function reduces to the Lorentz line shape and the resulting expression can be differentiated [Humlecek, 1979; ibid 1982].                f f f d Voigt 0 Lorentz 0 Doppler = (16)                                       f d L D Voigt 2 2 2 0 L 3/2 D L 0 = exp (17) This constitutes all of the terms which will appear in our expression. The cross-section temperate derivative can now formally be expressed using the chain rule                        dT df T dT f T dS T dT S T d T , , , = 1 1             (18) Meanwhile, the temperature derivatives of the line strength and the Lorentz width can be computed from their respective definitions.           S T  T c c T c T dT dQ T T Q T c E dT dS T                       2 2 0 2 0 2 0 2 2 1 exp / = 1 1 exp /    (19)   dT f d dT df f L L L L L L   1  = 1 2    1  (20) This expression, unfortunately, must be computed at the time that the cross-section database is created if we are to use the approach we have devised. This is because spectroscopic properties are used in the cross-section database creation stage but not carried through to the radiation transport (and retrieval) calculation. As a result, the size of the database we will need to carry has effectively doubled: for each cross-section we must now carry along its temperature derivative. This is not an onerous requirement as the original cross-section database is not excessively large. 4. OPTIMAL CHANNEL SELECTION Next we demonstrate how the new database aids in exploring concepts for our retrieval algorithm. Hyperspectral imagers acquire a substantial quantity of spectral data over a large frequency range and at a high resolution. The value of this data for retrievals varies spectrally. Here, we outline the procedure used to determine the optimal frequency ranges ("channels") to use for retrievals. Some data points are critically important to our analysis while others will provide little or no useful information. The use of a spectral filter (or mask) to separate the data into these two categories can both reduce computational time and increase the accuracy of the retrieval. This selection must be tailored to the target instrument. If we consider the ABB Bomem interferometer, our target instrument would be a 0.48217 wavenumber resolution over a relatively large spectral range. Many of the individual spectral frequencies (or channels) will not contain useful information. Perhaps those channels give radiation for molecules which are not of interest. Or, perhaps they are frequencies in which little or no atmospheric or HR-VOC radiation is expected. Finally, they may be channels in which the signal we are searching for is weak relative to other sources of radiation. Removing those channels with a mask function prior to applying a retrieval algorithm can be a way to improve the signal to noise ratio of our results. Further, the presence of absorption coefficient databases which span a much larger frequency space allows us test mask values for spectral ranges beyond the band center region. We define a "good" channel here as one in which measured radiance is dominated by a single molecular species. Such channels reduce the need to simultaneously retrieve parameters for multiple species. A good initial guess of the temperature and molecular concentration profiles can then be constructed by building up that guess one species at a time. As a final step, the retrieval algorithm will be run using the full data set with the composite single-species retrievals as initial conditions. This allows the algorithm to make small corrections for channels which molecules other than the dominant radiator make a small, but non-zero, contribution to the measured signal. Figure 1 prov retrieval algo instrument r designed to e to the compu vides a compa orithm on a l esolution, a v explore the opt utation, as is th Figu arison between logarithmic sc alue close to timal retrieval e temperature a ure 1. Selected our cross-sect ale. The selec but specificall approach for a and pressure ra d spectral chann tion calculation cted "instrume ly different th range of senso ange for which nels compared ns and the data ent" was artifi han the Bomem ors. Spectral re h to consider an to computed c a points select ficial. We cho m interferomet esolution is the nd the specific cross-sections. ted for inclusio se a 0.5 wave ter. The softw erefore simply selection crite on in our enumber ware was an input eria. We demonst calculations g designed for SHARC/SAM We illustrate obtain a tota temperature algorithms ar Figure 2. C We next exam comparison retrieval. Par We are unab series of pre using spectro the two code trate the retriev given by the S r the needs MM Atmosphe e the combined al radiance flu value (300 K re not expected omparison of t mine the tempe code's non-L rtition function le to compute e-computed cro oscopic databa s. Indeed, that F val algorithm b AMM code [D of this probl ere Generagor product of pat uctuation in tem K). The plots o d to adversely i the temperature FC erature derivat LTE treatment ns are stripped temperature de oss-sections. T ases which diff is what we see igure 3. Radian 5. DEM by applying it Dothe et al., 20 em [Panfili e (SAG) [Shroll th radiance and mperature due overlap and s impact our retr e derivative of CN subroutine ive of the optic means SAM d from spectros erivatives of lin These cross-sec fer from SAMM e in Figure 3. nce fluctuation MONSTRAT to determine 004] and the lin et al., 2011a] et al., 2003]. d the temperatu to the source how the same rieval. f the logarithm in the retrieva cal depth. This MM's computat scopic lines an nes at runtime ction derivativ M. We therefo ns in temperatu TION the densities a ne-by-line radi ]. The atmosp ure derivative o e term. This is e trend. Any of the source t al algorithm. s term is compu ational approac nd replaced wi as the required ves were comp ore expect clos ure due to the s and temperatur ation transport pheric profile of the logarithm s illustrated in subtle differen term computed uted at runtime ch diverges s ith their equiv d spectroscopy puted assumin se, but not exa source term. res from path t algorithm spe e is provided m of the source n Figure 2 for nces between d with SAMM e in SAMM. N substantially fr valent non-LTE y has been redu g LTE conditi act agreement radiance ecifically by the e term to a single the two and the Notice the from our E values. uced to a ions and between Finally, we i selected for layers are tw derivative fu temperature identical, agr Figure 4. logarithm We have out utilizes a we the radiation set with whi interaction. In the future, Spectral sele magnitude fo around a few Finally, the volumetric fl illustrate a sce optimal determ wo separate tem unction utilize (300 K) and a reement betwee Temperature d of the cross-se tlined the proce ell-established transport equa ich we perform , we see three ection of inter or our example w central frequ single line-of-lows of HR-VO enario consistin mination of wa mperatures, com d by the retri another layer c en the two algo derivatives of t ection (path-ind 6. CONC edure used to p fitting techniq ation within a l m the retrieval valuable impro rferometer dat e). Further per uencies. Next, -sight measure OCs through a ng of a three-ater content are mpare temperat ieval algorithm onsistent with orithms, as exp the cross-sectio dependent) and LUSIONS A predict molecu que coupled wi ine-by-line tre l has been out ovements to th ta reduces the rformance gain , the entire pro ments of insta defined volum -segment path e included. Th ture derivative m. The exam the temperatu pected. ons of water. T d the bottom ro section. AND FUTU ular concentrat ith analytic der atment. Finall tlined. This p his approach. T e computationa ns can be made ocess needs to antaneous CE me. from observer he derivative c es computed us mple shows re ure of a hot plu Top rows show ows show the te URE DIREC tions from obs rivatives obtai ly, a procedure procedure is fu The first is to i al load signif e by exploiting o be automated will need to b r to source. O computations, sing SAMM la esults for a la ume (800 K). w the temperatu emperature der CTION erved path rad ined from our e to reduce the ully automated improve comp ficantly (by m g the concentra d for turnkey u be used to con Only spectral c shown in Figu ayer calculation ayer near the We see close ure derivative o rivative of the diance. This pr directly differe overall size of d and requires putational perfo more than an ation of spectr usage by a non nstruct time-de channels ure 4 for ns to the ambient , but not of the cross-rocedure entiating f the data no user ormance. order of al points n-expert. ependent ACKNOWLEDGEMENTS Funding for this research comes from the Department of Energy through Contract No. DE-SC003373. Partial support for this research also comes from Spectral Sciences, Inc. internal research and development funding. Dr. Jeffrey A. Mercier of Sandia National Laboratory has provided substantial support to this project. We also acknowledge numerous substantive discussions of benefit to this effort with Mr. Russ Nettles of the Texas Commission on Environment Quality. We acknowledge our collaborators in this effort: Drs. Philip Smith and Jeremy Thornock (University of Utah), Dr. Kevin C. Gross and Mr. Michael R. Rhoby (Air Force Institute of Technology), and Dr. James Seebold (Chevron, ret.). REFERENCES Armstrong, B.H. Spectrum line profiles: the Voigt function. Journal of Quantitative Spectroscopy and Radiative Transfer, 7:61-88, (1967). Dothe, H., J.W. Duff, R. Panfili, JH. Gruninger, R.G. Kennett, P.K. Acharya, A. Berk, and L.S. Bernstein. Users' manual for SAMM2 (2004). Hartley, H.O. The modified gauss-newton method for the fitting of non-linear regression functions by least squares. Technometrics, 3(2):269-280, (1961). Humlicek, J. An efficient method for evaluation of the complex probability function: the voigt function and its derivatives. Journal of Quantitative Spectroscopy and Radiative Transfer, 21(4):309-313, (1979). Humlicek, J. Optimized computation of the voigt and complex probability functions. Journal of Quantitative Spectroscopy and Radiative Transfer, 27(4):437-444, (1982). Levenberg, K. A method for the solution of certain non-linear problems in least squares. Quart. Appl. Math., 2:164-168, (1944). Marquardt, D.W. An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial & Applied Mathematics, 11(2):431-441, (1963). More, Jorge J, Garbow, Burton S and Hillstrom, Kenneth E, User Guide for MINPACK, Argonne National Laboratory, report ANL-80-74 (1980). Panfili, R., Vujkovic-Cvijin, P, Tan, X, Kennett, R, Taylor, RS, Dothe, H, Bernstein, LS, Smith, P. Thornock, J, Gross, K and James Seebold, J, Remote Detection of Volatile Organic Compound Emissions from Combustion Flares, Proceedings of the Conference on Atmospheric Transmission Models 33 (2011). Panfili, R., H. Dothe, P. Vujkovic-Cvijin, X. Tan, R. Kennett, R. Taylor, and L. Bernstein, P. Smith, J. Thornock, K. Gross, J. Seebold, High-Fidelity Modeling of Flare Combustion and Emission Detection, American Flame Research Committee 2011 AFRC Combustion Symposium, Houston, Texas, September 18-21, (2011). Rothman, L.S., I.E. Gordon, A. Barbe, D. ChrisBenner, P.F. Bernath, M. Birk, V. Boudon, L.R. Brown, A. Campargue, J.-P. Champion, K. Chance, L.H. Coudert, V. Dana, V. M.Devi, S. Fally, J.-M. Flaud, R.R. Gamache, A. Goldman, D. Jacquemart, I. Kleiner, N. Lacome, W.J. Lafferty, J.-Y. Mandin, S.T. Massie, S.N. Mikhailenko, C.E. Miller, N. Moazzen-Ahmadi, O.V. Naumenko, A.V. Nikitin, J. Orphal, V.I. Perevalov, A. Perrin, A. Predoi-Cross, C.P. Rinsland, M. Rotger, M.Simeckova, M.A.H. Smith, K. Sung, S.A. Tashkun, J. Tennyson, R.A. Toth, A.C. Vandaele, and J. VanderAuwera, "The HITRAN 2008 molecular spectroscopic database," Journal of Quanitative Spectrosscopy and Radiative.Transfer 110, 533-572 (2009). Rothman, L.S., I.E. Gordon, R.J. Barber, H. Dothe, R.R. Gamache, A. Goldman, V.I. Perevalov, S.A. Tashkun, and J. Tennyson. Hitemp, the high-temperature molecular spectroscopic database. Journal of Quantitative Spectroscopy and Radiative Transfer, 111(15):2139-2150, (2010). Seebold, J., P. Gogoloek, J. Pohl, and R. Schwartz. Practical implications of prior research on today's outstanding flare emissions questions and a research program to answer them. In AFRC-JFRC 2004 Joint International Combustion Symposium, Maui, HI, (2004). Seebold, J. The origin and fate of toxic combustion byproducts in refinery heaters: Research to enable efficient compliance with the clean air act. Technical Report, Petroleum Environmental Research Forum Project 92-19, Final Report, 1997. Seebold, J.G. and R.T. Waibel. Products of incomplete combustion (pic) from petroleum, petrochemical & chemical sector process heaters and industrial boilers. In 10th International Congress on Combustion By-Products and Their Health Effects, Ischia, Italy, (June 17-20, 2007), 2007. Shroll, R.M, S. Adler-Golden, J.W. Duff, and J.H. Brown. User's manual for SAG-2 SHARC/SAMM atmosphere generator (2003).