A polarized discrete ordinatescattering modelfor radiative transfer simulationsin spherical atmosphereswith thermal source

Vom Fachbereich für Physik und Elektrotechnikder Universität Bremenzur Erlangung des akademischen Grades einesDoktor der Naturwissenschaften (Dr. rer. nat.)genehmigte Dissertation

vonDipl. Phys. Claudia Emde

Dezember 2004

Berichte aus dem Institut für Umweltphysik – Band 25herausgegeben von:

Dr. Georg HeygsterUniversität Bremen, FB 1, Institut für Umweltphysik,Postfach 33 04 40, D-28334 BremenURL http://www.iup.physik.uni-bremen.deE-Mail iupsekr@uni-bremen.deDie vorliegende Arbeit ist die inhaltlich unveränderte Fassung einer Dissertation,die im Dezember 2004 dem Fachbereich Physik/Elektrotechnik der UniversitätBremen vorgelegt und von Prof. Dr. Klaus Künzi sowie Prof. Dr. Clemens Simmerbegutachtet wurde. Das Promotionskolloquium fand am 28. Januar 2005 statt.

Bibliografische Information Der Deutschen BibliothekDie Deutsche Bibliothek verzeichnet diese Publikation in derDeutschen Nationalbibliografie; detaillierte bibliografische Datensind im Internet über http://dnb.ddb.de abrufbar.

c© Copyright 2005 Logos Verlag Berlin

Alle Rechte vorbehalten.

ISBN 3-8325-0855-4 ISSN 1615-6862

Logos Verlag BerlinComeniushofGubener Straße 47D-10243 BerlinTelefon (0 30) 42 85 10 90URL http://www.logos-verlag.de

Layout: Lothar Meyer-Lerbs, Bremen

Contents

Publications 10

Preface 15

1 Theoretical background 191.1 Basic definitions 191.2 The Stokes parameters 201.3 Single particle scattering 23

1.3.1 Definition of the amplitude matrix 241.3.2 Phase matrix 251.3.3 Extinction matrix 261.3.4 Absorption vector 261.3.5 Optical cross sections 28

1.4 Particle Ensembles 291.4.1 Single scattering approximation 29

1.5 Radiative transfer equation 32

2 ARTS 372.1 History 372.2 Definition of the atmosphere 38

2.2.1 Atmospheric dimensionality 382.2.2 The cloud box 39

2.3 Radiative transfer calculations 402.3.1 Propagation paths 412.3.2 Radiative background 422.3.3 Clear sky radiative transfer 43

2.4 Scattering 442.5 Gas absorption 442.6 Definition of clouds and atmospheric fields 44

3

4 Contents

2.7 Unit conventions 45

3 Clouds as scattering media 473.1 Microphysics of clouds 473.2 Coordinate systems 503.3 Computation methods 52

3.3.1 Rayleigh scattering 523.3.2 Lorentz-Mie theory for scattering by spherical

particles 543.3.3 T-matrix method 563.3.4 Further methods 63

3.4 Single scattering properties ARTS 653.5 Particle size distributions 69

3.5.1 Mono-disperse particle distribution 693.5.2 Gamma size distribution 703.5.3 McFarquhar and Heymsfield parametrization 71

4 The DOIT scattering model 734.1 The discrete ordinate iterative method 73

4.1.1 Radiation field 734.1.2 Vector radiative transfer equation solution 744.1.3 Scalar radiative transfer equation solution 784.1.4 Single scattering approximation 79

4.2 Sequential update 804.2.1 Up-looking directions 814.2.2 Down-looking directions 824.2.3 Limb directions 83

4.3 Grid optimization and interpolation methods 844.3.1 Zenith angle grid optimization 854.3.2 Interpolation methods 874.3.3 Error estimates 88

5 Scattering model intercomparisons 915.1 FM2D - A 2D pseudo-spherical model 91

5.1.1 The pseudo-spherical approach 915.1.2 Clear sky comparison 935.1.3 Comparison for cloudy scenarios 94

Contents 5

5.2 KOPRA - A single scattering model for the IR 985.2.1 Zero- and single scattering solutions 985.2.2 Definition of scenarios 1005.2.3 Results for case ω0 = 0.24 1025.2.4 Results for case ω0 = 0.84 1045.2.5 Summary and discussion 105

5.3 ARTS - Monte Carlo Scattering Model 1085.3.1 The Monte Carlo approach 1085.3.2 Setup 1105.3.3 Results 1125.3.4 Discussion 1145.3.5 Summary and conclusions 118

6 Microwave limb spectra 1196.1 General setup for the simulations 119

6.1.1 Atmosphere 1206.1.2 Sensor setup 1216.1.3 Numerical setup 121

6.2 Definition of cloud scenarios 1226.3 Results 123

6.3.1 Impact of particle size 1236.3.2 Impact of cloud altitude 1256.3.3 Impact of particle shape 1266.3.4 Correlated IMC and Reff 1286.3.5 Spectra for different frequency bands 1286.3.6 Comparison with nadir radiances 131

6.4 Discussion 136

7 Simulation of polarized radiances 1397.1 Model simulations in a 1D spherical atmosphere 139

7.1.1 Scattering and polarization signal for differentparticle sizes 140

7.1.2 Effect of particle shape 1447.1.3 Scalar simulations 144

7.2 3D box type cloud model simulations 1487.2.1 Performance 151

7.3 Conclusions 153

6 Contents

8 Thin layer cirrus study for EOS MLS 1578.1 Setup 157

8.1.1 Selection of frequencies 1578.1.2 Particle size distribution function 1588.1.3 Included species for absorption coefficient calculations 1588.1.4 Definition of realistic cloud parameters 159

8.2 Clear sky and cloudy radiances 1598.3 Comparison between 1D and 3D simulations 1608.4 Sensitivity study 163

8.4.1 Dependence on ice mass content 1638.4.2 Dependence on aspect ratio 164

8.5 Combination of polarized channels 1678.6 Conclusions and outlook 168

9 Summary, conclusions and outlook 171

Appendix 175

Appendix 177

A Literature review 177

B Derivations 181B.1 Solution of approximated VRTE 181B.2 Coordinate system transformation 182

B.2.1 Phase matrix transformation 182B.2.2 Extinction matrix and absorption vector 184

C List of acronyms 189

D List of symbols 191

Abstract

This work describes the development of the new discrete ordinatescattering algorithm, which is a part of the Atmospheric RadiativeTransfer Simulator (ARTS). Furthermore, applications of the algo-rithm, which was implemented to study for example the influence ofcirrus clouds on microwave limb sounding, are presented.

The model development requires as a theoretical basis the electro-magnetic scattering theory. The basic quantities are defined and differ-ent methods to compute single scattering properties of small particlesare discussed. The phenomenological derivation of the vector radiativetransfer equation, which is the basic equation of the model, is outlined.In order to represent clouds as scattering media in radiative transfermodels, information about their micro-physical state is required asan input for calculating the scattering properties. The micro-physicalstate of a cloud is defined by the phase of the cloud particles, theparticle size and shape distributions, the particle orientation, the icemass or the liquid water content, and the temperature.

The model uses the Discrete Ordinate ITerative (DOIT) method tosolve the vector radiative transfer equation. The implementation of adiscrete ordinate method is challenging due to the spherical geome-try of the model atmosphere, which is required for the simulation oflimb radiances. The involved numerical issues, grid optimization andinterpolation methods, are discussed.

The new scattering algorithm was compared to three other models,which were developed during the same time period as the DOIT al-gorithm. Overall, the agreement between the models was very good,giving confidence in new models.

Scattering simulations are presented for limb- and down-looking ge-ometries, for one-dimensional and three-dimensional spherical atmo-

7

8

spheres. They were performed for the frequency bands of the Mil-limeter Wave Acquisitions for Stratosphere/Troposphere ExchangeResearch (MASTER) instrument, and for selected frequencies of theEarth Observing System Microwave Limb Sounder (EOS MLS). Thesimulations show the impact of cloud particle size, shape and orien-tation on the brightness temperatures and on the polarization of mi-crowave radiation in the atmosphere. The cloud effect is much largerfor limb radiances than for nadir radiances. Particle size is a veryimportant parameter in all of the simulations. The polarization sig-nal is small for simulations with randomly oriented particles whereasfor horizontally aligned particles with random azimuthal orientationthe polarization signal is significant. Moreover, the effect of particleshape is only relevant for oriented cloud particles. The simulationsshow that it is essential to use a three-dimensional scattering modelfor inhomogeneous cloud layers.

Publications

The work described in this text has given rise to a number of publi-cations:

Journal Articles1. A detailed description of the DOIT scattering model including sev-

eral examples is published in:Emde, C., S. A. Buehler, C. Davis, P. Eriksson, Sreerekha T. R.and C. Teichmann, A polarized discrete ordinate scattering modelfor simulations of limb and nadir longwave measurements in 1D/3Dspherical atmospheres, J. Geophys. Res., in press 2004.

2. First limb scattering simulations for MASTER frequency bands arepresented in:Emde, C., S. A. Buehler, P. Eriksson and T. R. Sreerekha (2004),The effect of cirrus clouds on microwave limb radiances, J. Atmos.Res., 72(1–4), 383–401, doi:10.1016/j.atmosres.2004.03.023.

3. The intercomparison of the ARTS-DOIT model with the radiativetransfer model KOPRA is described in:Höpfner, M. and C. Emde (2005), Comparison of single and multiplescattering approaches for the simulation of limb-emission observa-tions in the mid-IR, J. Quant. Spectrosc. Radiat. Transfer, 91(3),275–285, doi:10.1016/j.jqsrt.2004.05.066.

4. The following article includes a description of the Monte Carlo scat-tering model, which has been implemented in ARTS by C. Davis,and several calculations using this model:Davis, C., C. Emde and R. Harwood, A 3D polarized reversed Monte

9

10 Publications

Carlo radiative transfer model for mm and sub-mm passive remotesensing in cloudy atmospheres, IEEE T. Geosci. Remote, in press,2004.

5. A prototype of the DOIT scattering module and a sensitivity studyof cloud parameters on nadir radiances is presented in:

Sreerekha, T. R., S. A. Buehler and C. Emde (2002), A simple newradiative transfer model for simulating the effect of cirrus clouds inthe microwave spectral region, J. Quant. Spectrosc. Radiat. Trans-fer, 75, 611–624.

6. A paper including an explanation of polarization induced by spher-ical particles using the DOIT model has been submitted:

Teichmann, C., S. A. Buehler and C. Emde Understanding the po-larization signal of spherical particles for microwave limb radiancesJ. Geophys. Res., submitted 2004.

7. The ARTS clear sky model was intercompared with several radiativetransfer models. The results are published in the following paper:

Melsheimer, C., C. Verdes, S. A. Buehler, C. Emde, P. Eriksson,D. G. Feist, S. Ichizawa, V. O. John, Y. Kasai, G. Kopp, N. Koulev,T. Kuhn, O. Lemke, S. Ochiai, F. Schreier, T. R. Sreerekha, M. Su-zuki, C. Takahashi, S. Tsujimaru and J. Urban, Intercomparison ofGeneral Purpose Clear Sky Atmospheric Radiative Transfer Mod-els for the Millimeter/Submillimeter Spectral Range, Radio Sci., inpress, 2004.

Technical Reports8. The ARTS user guide also includes a detailed theoretical description

and furthermore information about the usage of the model:

Eriksson, P., S. A. Buehler, C. Emde, T. R. Sreerekha, C. Melsheimerand O. Lemke (2004), ARTS-1-1 User Guide, University of Bre-men, 308 pages, regularly updated versions available at www.sat.uni-bremen.de/arts/.

Publications 11

9. A large part of the work was carried out in the context of ESCTECstudies: The UTLS study consisted of two tasks. The topic of task1 was “2D retrievals of cloud-free scenes”. Task 2 was about 2D re-trievals in the presence of clouds, for this part the DOIT model wasdeveloped and applied for simulations as well as for model compar-isons. The executive summary includes the most important findingsand conclusions:

Kerridge, B., V. Jay, J. Reburn, R. Siddans, B. Latter, F. Lama,A. Dudhia, D. Grainger, A. Burgess, M. Höpfner, T. Steck, C. Emde,P. Eriksson, M. Ekström, A. Baran and M. Wickett (2004), Con-sideration of mission studying chemistry of the UTLS, ExecutiveSummary, ESTEC Contract No 15457/01/NL/MM.

10. The results of task 2 are compiled in:

Kerridge, B., V. Jay, J. Reburn, R. Siddans, B. Latter, F. Lama,A. Dudhia, D. Grainger, A. Burgess, M. Höpfner, T. Steck, G. Stiller,S. Buehler, C. Emde, P. Eriksson, M. Ekström, A. Baran andM. Wickett (2004), Consideration of mission studying chemistry ofthe UTLS, Task 2 Report, ESTEC Contract No 15457/01/NL/MM.

11. The final report includes the major results of task 1 and task 2:

Kerridge, B., V. Jay, J. Reburn, R. Siddans, B. Latter, F. Lama,A. Dudhia, D. Grainger, A. Burgess, M. Höpfner, T. Steck, G. Stiller,S. Buehler, C. Emde, P. Eriksson, M. Ekström, A. Baran andM. Wickett (2004), Consideration of mission studying chemistry ofthe UTLS, Final Report, ESTEC Contract No 15457/01/NL/MM.

12. For the RT study a literature review of radiative transfer modelsfor the microwave region which include scattering was performed:

Claudia Emde and Sreerekha T. R. (2004), Development of a RTmodel for frequencies between 200 and 1000 GHz, WP1.2 ModelReview, ESTEC Contract No AO/1-4320/03/NL/FF.

12 Publications

Articles in Conference Proceedings13. Preliminary results obtained by using the DOIT scattering module

are shown in:Emde, C., S. Buehler, Sreerekha T. R. (2003), Modeling polarizedmicrowave radiation in a 3D spherical cloudy atmosphere, In: Elec-tromagnetic and Light Scattering – Theory and applications VII,Edited by Wriedt, T., Universität Bremen, ISBN 3-88722-579-1.

14. Nadir simulations to estimate the cloud effect for AMSU-B channelsare published in:Sreerekha, T. R., C. Emde and S. A. Buehler, Using a new Ra-diative Transfer Model to estimate the Effect of Cirrus Clouds onAMSU-B Radiances, In: Twelfth International TOVS Study Con-ference (ITSC – XII), Lorne, Australia, February 2002.

15. The DOIT model was compared to the RTTOVSCAT model andto AMSU radiances. The results are presented in the conferenceproceedings:English, S. J., U. O’Keeffe, T. R. Sreerekha, S. A. Buehler, C. Emdeand A. M. Doherty (2003), A Comparison of RTTOVSCATT withARTS and AMSU Observations, Using the Met Office MesoscaleModel Short Range Forecasts of Cloud Ice and Liquid Water, In:Thirteenth International TOVS Study Conference (ITSC – XIII),St. Adele, Montreal, Canada.

Preface

To improve existing climate models it is very important to extend theknowledge about cirrus cloud parameters, as such clouds cover morethan 20% of the globe (Wang et al., 1996) and play an importantrole in the Earth’s radiation budget (Arking, 1991). Depending oncloud altitude and micro-physical properties, clouds can either causewarming or cooling at the Earth’s surface. So far clouds are not welltreated in Global Climate Models (GCM) because of uncertaintiesconcerning the properties of cirrus clouds and because of the complexinteraction between radiation, micro-physics and dynamics in theseclouds. Moreover it is essential to consider clouds for the evaluation oflimb measurements of trace gases in the upper troposphere. Clouds,especially cirrus, with particle sizes exceeding microwave wavelengths,can severely disturb trace gas measurements. On the other hand it ispossible to obtain cloud information from microwave limb radiancesaffected by cirrus clouds. This requires a radiative transfer model thatcan simulate the scattering effect of cirrus clouds.

In particular the effective radius Reff of cloud particles is impor-tant for the radiative properties of clouds. For a given frequency, Reff

largely determines the relation between ice mass content (IMC) andcloud optical thickness (Evans et al., 1998). Parameterizations of Reff

have been retrieved from combined lidar and radar reflectivity (Dono-van, 2003) or from observations made in situ using aircraft mountedinstruments (e.g., Kinne et al., 1997; McFarquhar and Heymsfield,1997). The Submillimeter-Wave Cloud Ice Radiometer (SWCIR) tofly on an aircraft has been developed to retrieve upper troposphericIMC and Reff (Evans et al., 2002).

Satellite remote sensing techniques in the thermal infrared can onlybe applied for thin cirrus clouds consisting of small ice particles as

13

14 Preface

saturation is reached for moderate optical depths (Stubenrauch et al.,1999). Only ice particle properties of the uppermost cloud layers canbe measured. Disadvantages of visible and near-infrared solar reflec-tion methods include that they cannot measure low optical depthclouds over brighter land surfaces.

There are several studies about the sensitivity of cirrus cloudson microwave nadir radiances (e.g., Evans et al., 1998; Skofronik-Jackson et al., 2002). They show that the brightness temperaturedepression depends strongly on particle size and IMC. Passive nadir-viewing techniques cannot sufficiently resolve the vertical distributionof IMC. Millimeter-wave limb sounding is a well established techniquefor the observation of atmospheric trace gases in the stratosphereand upper troposphere. This technique can provide higher resolu-tions than nadir techniques for the same frequencies. Instruments us-ing this technique are the Earth Observing System Microwave LimbSounder (EOS MLS) (Waters et al., 1999), the Millimeter Atmo-spheric Sounder (MAS) (Hartmann et al., 1996) and the MillimeterWave Acquisitions for Stratosphere/Troposphere Exchange Research(MASTER) instrument (Buehler, 1999). Recently, instruments havemoved towards higher frequencies into the submillimeter-wave region,examples of this type of instrument are Odin-SMR (Murtagh et al.,2002) and the Superconduction Submillimeter-Wave Limb EmissionSounder (SMILES) (Buehler et al., 2005b).

A number of well established radiative transfer models exist forthe clear sky case, notably the public domain Atmospheric RadiativeTransfer Simulator (ARTS) (Buehler et al., 2005a), which was takenas the platform for the new scattering model described in this the-sis. The model development is a challenging task for various reasons:Firstly, cloud coverage is vertically and horizontally inhomogeneouswhich implies that a three-dimensional (3D) model is unavoidable forthe simulation of realistic cases. Especially for limb measurements,the 3D spherical geometry is required as the observed region in theatmosphere has a horizontally large extent. However, for largely ex-tended thin cirrus clouds, it makes sense to use a one-dimensional(1D) model, because it can be much more efficient compared to a full3D model. Secondly, cirrus clouds consist of particles of different sizes

Preface 15

and shapes. Since particle scattering due to non-spherical particlesleads to polarization effects (Czekala and Simmer, 1998), the vectorradiative transfer equation (VRTE) has to be used in the model toobtain the full Stokes vector, not just the intensity of the radiation.

Liquid water clouds are not so problematic, because liquid waterdrops mainly act as absorbers, not as scatterers. Cirrus clouds, onthe other hand, have a low absorption coefficient (see for exampleMishchenko et al., 2002) and a rather large scattering coefficient.Aerosol scattering needs to be considered in the infra-red. Molecu-lar Rayleigh scattering, though important for optical wavelength, canbe neglected at microwave and infra-red wavelengths.

A survey of existing freely available radiative transfer models yieldednone that were well-suited to the requirements described above. Forinstance the 3D Monte Carlo models described in Liu et al. (1996) andRoberti et al. (1994) are only applicable for 3D-cartesian atmospheres.The 3D discrete ordinate models SHDOM (Evans, 1998) and VDOM(Haferman et al., 1997) also assume a cartesian geometry. For thisreason they are not applicable for limb simulations. Other discrete or-dinate models, for example MWMOD (Simmer, 1993) and VDISORT(Schulz and Stamnes, 2000), use one-dimensional (1D) plane-parallelgeometries. Another well known method is the Eddington approxima-tion (e.g., Kummerow, 1993), which is also not well-suited to the limbsounding problem, as it is only valid in plane-parallel atmospheres. Asimple 1D plane-parallel model using a prototype of the iterative so-lution method described in this thesis is presented in Sreerekha et al.(2002).

In the new version of ARTS two scattering methods have been im-plemented: a backward Monte Carlo Method (Davis et al., 2004) andthe DOIT (Discrete Ordinate ITerative) (Emde et al., 2004a) methodbeing presented in this work. Both methods work in 3D spherical at-mospheres and both can simulate polarization effects due to asphericalparticles. The DOIT method works also in 1D spherical atmospheres.The implementation of the DOIT method is very similar to discreteordinate method (DOM) implementations for instance in SHDOM orVDOM. The originality of the DOIT method is, that the DOM hasbeen adapted to a spherical geometry, which is essential for the simu-

16 Preface

lation of limb radiances. The model can be applied in the microwaveand in the infrared wavelength regions.

The present thesis consists of nine chapters. Chapter 1 gives an in-troduction to the theoretical concepts of the radiative transfer theoryfor scattering media. Basic quantities are defined and an outline of thephenomenological derivation of the vector radiative transfer equationis given. Chapter 2 introduces concepts and definitions of ARTS, whichwas used as a platform to implement the DOIT algorithm. Chapter 3gives a brief overview of cloud micro-physics and it introduces differ-ent methods to calculate scattering properties of cloud particles. It isshown that the T-matrix method is the most appropriate to be usedfor the new scattering model.

The DOIT algorithm is described in detail in Chapter 4, whichstarts with the theoretical basis of discrete ordinates and afterwardsexplains the numerical optimizations, which were necessary for effi-ciency reasons. In Chapter 5 the 1D comparisons with the modelsFM2D and KOPRA are shown. Here ARTS-DOIT was used as a ref-erence model, since it is the more general and more accurate model.Furthermore a 3D comparison with the ARTS Monte Carlo model ispresented. Note that the two scattering models presented here are thefirst models which are able to simulate polarization in a 3D sphericalatmosphere in the microwave region.

Chapter 6 presents 1D simulations for the MASTER instrument,where the effect of cloud parameters like effective radius and ice massconstant is investigated. First simulations using the full capabilitiesof the new model, i.e., polarization and 3D geometry, are shown inChapter 7. In Chapter 8 simulations for the EOS MLS instrument arepresented. Scattering and polarization of thermal radiation in a thinlayer tropical cirrus cloud are investigated.

The final Chapter 9 consists of the overall summary and of conclu-sions.

1 Theoretical background

This chapter introduces the theoretical background which is essen-tial to develop a radiative transfer model including scattering. Thetheory is based on concepts of electrodynamics, starting from theMaxwell equations. An elementary book for electrodynamics is writ-ten by Jackson (1998). For optics and scattering of radiation by smallparticles the reader may refer for instance to van de Hulst (1957)and Bohren and Huffman (1998). The notation used in this chapteris mostly adapted from the book “Scattering, Absorption, and Emis-sion of Light by Small Particles” by Mishchenko et al. (2002). Severallengthy derivations of formulas, which are not shown in detail here,can also be found in this book. The purpose of this chapter is to pro-vide definitions and give ideas, how these definitions can be derivedusing principles of electromagnetic theory. For the derivation of theradiative transfer equation an outline of the traditional phenomeno-logical approach is given.

1.1 Basic definitionsFrom the Maxwell equations one can derive the formula for the electro-magnetic field vector E of a plane electromagnetic wave propagatingin a homogeneous medium without sources:

E(r, t) = E0 exp(−ω

cmIn · r

)exp

(iω

cmRn · r − iωt

), (1.1)

where E0 is the amplitude of the electromagnetic wave in vacuum,c is the speed of light in vacuum, ω is the angular frequency, r isthe position vector and n is a real unit vector in the direction ofpropagation. The complex refractive index m is

m = mR + imI = c√

εµ, (1.2)

17

18 1 Theoretical background

where mR is the non-negative real part and mI is the non-negativeimaginary part. Furthermore µ is the permeability of the medium andε the permittivity. For a vacuum, m = mR = 1. The imaginary partof the refractive index, if it is non-zero, determines the decay of theamplitude of the wave as it propagates through the medium, which isthus absorbing. The real part determines the phase velocity v = c/mR.The time-averaged Poynting vector P (r), which describes the flow ofelectromagnetic energy, is defined as

P (r) =12Re(〈E(r)〉 × 〈H∗(r)〉

), (1.3)

where H is the magnetic field vector and the ∗ denotes the complexconjugate. The Poynting vector for a homogeneous wave is given by

〈P (r)〉 =12Re(√

ε

µ

)|E0|2 exp

(−2

ω

cmIn · r

)n. (1.4)

Equation (1.4) shows that the energy flows in the direction of propa-gation and its absolute value I(r) = |〈P (r)〉|, which is usually calledintensity (or irradiance), is exponentially attenuated. Rewriting Equa-tion (1.4) gives

I(r) = I0 exp(−αpn · r), (1.5)

where I0 is the intensity for r = 0. The absorption coefficient αp is

αp = 2ω

cmI =

4πmI

λ=

4πmIν

c, (1.6)

where λ is the free-space wavelength and ν the frequency. Intensityhas the dimension of monochromatic flux [energy/(area × time)].

1.2 Definition of the Stokes parametersSensors usually do not measure directly the electric and the magneticfields associated with a beam of radiation. They measure quantitiesthat are time averages of real-valued linear combinations of productsof field vector components and have the dimension of intensity. Ex-amples of such observable quantities are the Stokes parameters. Fig-ure 1.1 shows the coordinate system used to describe the direction of

1.2 The Stokes parameters 19

propagation n and the polarization state of a plane electromagneticwave.

z

x

y

n

φ

θ

φ

θ

O

Figure 1.1: Coordinate system to describe the direction of propagationand the polarization state of a plane electromagnetic wave (adapted fromMishchenko).

The unit vector n can equivalently be described by a couplet (θ, φ),where θ ∈ [0, π] is the polar (zenith) angle and φ ∈ [0, 2π) is theazimuth angle. The electric field at the observation point is given byE = Eθ + Eφ, where Eθ and Eφ are the θ- and φ-components of theelectric field vector. Eθ lies in the meridional plane, which is the planethrough n and the z-axis, and Eφ is perpendicular to this plane. Eθ

and Eφ are often called Ev and Eh in the microwave remote sensingliterature.

20 1 Theoretical background

The Stokes parameters are defined as follows:

I = 12

√εµ (EvE∗

v + EhE∗h), (1.7)

Q = 12

√εµ (EvE∗

v − EhE∗h), (1.8)

U = − 12

√εµ (EvE∗

h + EhE∗v), (1.9)

V = i 12

√εµ (EhE∗

v − EvE∗h). (1.10)

They are commonly defined as a 4×1 column vector I, which is knownas the Stokes vector. Since the Stokes parameters are real-valued andhave the dimension of intensity, they can be measured directly withsuitable instruments. The Stokes parameters are a complete set ofquantities needed to characterize a plane electromagnetic wave. Theycarry information of the complex amplitudes and the phase difference.The first Stokes parameter I is the intensity and the other componentsQ, U and V describe the polarization state of the wave. The Stokes pa-rameters of a plane monochromatic wave are related by the quadraticidentity

I2 = Q2 + U2 + V 2. (1.11)

The definition of a monochromatic plane wave implies that the com-plex amplitude E0 and the phase differences are constant. In the caseof natural radiation the amplitudes and phases fluctuate, since theradiation originates from several sources that do not emit radiationcoherently, and since the emission from one source usually has veryshort coherence times. This means that we usually have a superposi-tion of radiation from several incoherent sources, and that the polar-ization state of the radiation from each source fluctuates as well. Suchfluctuations have time scales that are longer than the period (2π/ω)of the oscillation, but that are still shorter than the integration timeof the instrument that measures the radiation. Thus, the instrumentmeasures an incoherent superposition of time averages over the fluc-tuating polarization. If the fluctuations are not completely random,the radiation is called partially polarized.

Since the different sources and/or emission events are assumed to

1.3 Single particle scattering 21

be incoherent, the Stokes parameters can simply be added up:

I =∑

i

Ii, Q =∑

i

Qi, U =∑

i

Ui, V =∑

i

Vi. (1.12)

The equality Equation (1.11) still holds for each contribution i, butfor the resulting I, Q, U , V , we have in general the inequality

I2 ≥ Q2 + U2 + V 2. (1.13)

The degree of polarization p is defined as

p =

√Q2 + U2 + V 2

I. (1.14)

For completely polarized radiation, Q2 + U2 + V 2 = I2, thus p = 1,and for unpolarized radiation, Q = U = V = 0, thus p = 0.

In addition to the degree of polarization, p, the degree of linearpolarization is defined as

plin =

√Q2 + U2

I, (1.15)

and the the degree of circular polarization is defined as

pcirc =V

I. (1.16)

1.3 Scattering, absorption and thermalemission by a single particle

A parallel monochromatic beam of electromagnetic radiation propa-gates in vacuum without any change in its intensity or polarizationstate. A small particle, which is interposed into the beam, can causeseveral effects:Absorption: The particle converts some of the energy contained in

the beam into other forms of energy.Elastic scattering: Part of the incident energy is extracted from the

beam and scattered into all spatial directions at the frequency ofthe incident beam. Scattering can change the polarization state ofthe radiation.

22 1 Theoretical background

Extinction: The energy of the incident beam is reduced by an amountequal to the sum of absorption and scattering.

Dichroism: The change of the polarization state of the beam as itpasses a particle.

Thermal emission: If the temperature of the particle is non-zero,the particle emits radiation in all directions over a large frequencyrange.

The beam is an oscillating plane magnetic wave, whereas the parti-cle can be described as an aggregation of a large number of discreteelementary electric charges. The incident wave excites the charges tooscillate with the same frequency and thereby radiate secondary elec-tromagnetic waves. The superposition of these waves gives the totalelastically scattered field.

One can also describe the particle as an object with a refractiveindex different from that of the surrounding medium. The presence ofsuch an object changes the electromagnetic field that would otherwiseexist in an unbounded homogeneous space. The difference of the totalfield in the presence of the object can be thought of as the field scat-tered by the object. The angular distribution and the polarization ofthe scattered field depend on the characteristics of the incident fieldas well as on the properties of the object as its size relative to thewavelength and its shape, composition and orientation.

1.3.1 Definition of the amplitude matrix

For the derivation of a relation between the incident and the scatteredelectric field we consider a finite scattering object in the form of asingle body or a fixed aggregate embedded in an infinite homogeneous,isotropic and non-absorbing medium. We assume that the individualbodies forming the scattering object are sufficiently large that they canbe characterized by optical constants appropriate to bulk matter, notto optical constants appropriate for single atoms or molecules. Solvingthe Maxwell equations for the internal volume, which is the interiorof the scattering object, and the external volume one can derive aformula, which expresses the total electric field everywhere in space

1.3 Single particle scattering 23

in terms of the incident field and the field inside the scattering object.Applying the far field approximation gives a relation between incidentand scattered field, which is that of a spherical wave. The amplitudematrix S(nsca, ninc) includes this relation:(

Escaθ (rnsca)

Escaφ (rnsca)

)=

eikr

rS(nsca, ninc)

(Einc

0θ

Einc0φ

). (1.17)

The amplitude matrix depends on the directions of incident ninc andscattering nsca as well as on size, morphology, composition, and orien-tation of the scattering object with respect to the coordinate system.The distance between the origin and the observation point is denotedby r and the wave number of the external volume is denoted by k.

The amplitude matrix provides a complete description of the scat-tering pattern in the far field zone. The amplitude matrix explicitlydepends on φinc and φsca even when θinc and/or θsca equal 0 or π.

1.3.2 Phase matrix

The phase matrix Z describes the transformation of the Stokes vectorof the incident wave into that of the scattered wave for scatteringdirections away from the incidence direction (nsca 6= ninc),

Isca(rnsca) =1r2

Z(nsca, ninc)I inc. (1.18)

The 4×4 phase matrix can be written in terms of the amplitude matrixelements for single particles (Mishchenko et al., 2002). All elementsof the phase matrix have the dimension of area and are real. As theamplitude matrix, the phase matrix depends on φinc and φsca evenwhen θinc and/or θsca equal 0 or π. In general, all 16 elements of thephase matrix are non-zero, but they can be expressed in terms of onlyseven independent real numbers. Four elements result from the moduli|Sij | (i, j = 1, 2) and three from the phase-differences between Sij .If the incident beam is unpolarized, i.e., I inc = (I inc, 0, 0, 0)T , thescattered light generally has at least one non-zero Stokes parameter

24 1 Theoretical background

other than intensity:

Isca = Z11Iinc, (1.19)

Qsca = Z21Iinc, (1.20)

U sca = Z31Iinc, (1.21)

V sca = Z41Iinc. (1.22)

This is the phenomena is traditionally called “polarization”. The non-zero degree of polarization Equation (1.14) can be written in terms ofthe phase matrix elements

p =

√Z2

21 + Z231 + Z2

41

Z11. (1.23)

1.3.3 Extinction matrix

In the special case of the exact forward direction (nsca = ninc) theattenuation of the incoming radiation is described by the extinctionmatrix K. In terms of the Stokes vector we get

I(rninc)∆S = I inc∆S −K(ninc)I inc + O(r−2). (1.24)

Here ∆S is a surface element normal to ninc. The extinction matrixcan also be expressed explicitly in terms of the amplitude matrix. Ithas only seven independent elements. Again the elements depend onφinc and φsca even when the incident wave propagates along the z-axis.

1.3.4 Absorption vector

The particle also emits radiation if its temperature T is above zeroKelvin. According to Kirchhoff’s law of radiation the emissivity equalsthe absorptivity of a medium under thermodynamic equilibrium. Theenergetic and polarization characteristics of the emitted radiationare described by a four-component Stokes emission column vectora(r, T, ω). The emission vector is defined in such a way that the netrate, at which the emitted energy crosses a surface element ∆S normal

1.3 Single particle scattering 25

to r at distance r from the particle at frequencies from ω to ω + ∆ω,is

W e =1r2

a(r, T, ω)B(T, ω)∆S∆ω, (1.25)

where W e is the power of the emitted radiation and B is the Planckfunction. In order to calculate a we assume that the particle is placedinside an opaque cavity of dimensions large compared to the par-ticle and any wavelengths under consideration. We have thermody-namic equilibrium if the cavity and the particle are maintained at theconstant temperature T . The emitted radiation inside the cavity isisotropic, homogeneous, and unpolarized. We can represent this radi-ation as a collection of quasi-monochromatic, unpolarized, incoherentbeams propagating in all directions characterized by the Planck black-body radiation

B(T, ω)∆S∆Ω =~ω3

2π2c2[exp

(~ω

kBT

)− 1]∆S ∆Ω, (1.26)

where ∆Ω is a small solid angle about any direction, ~ is the Planckconstant divided by 2π, and kB is the Boltzmann constant. The black-body Stokes vector is

Ib(T, ω) =

B(T, ω)

000

. (1.27)

For the Stokes emission vector, which we also call particle absorptionvector, we can derive

api (r, T, ω) = Ki1(r, ω)−

∫4π

dr′Zi1(r, r′, ω), i = 1, . . . , 4.

(1.28)

This relation is a property of the particle only, and it is valid for anyparticle, in thermodynamic equilibrium or non-equilibrium.

26 1 Theoretical background

1.3.5 Optical cross sections

The optical cross-sections are defined as follows: The product of thescattering cross section Csca and the incident monochromatic energyflux gives the total monochromatic power removed from the incidentwave as a result of scattering into all directions. The product of theabsorption cross section Cabs and the incident monochromatic energyflux gives the power which is removed from the incident wave by ab-sorption. The extinction cross section Cext is the sum of scattering andabsorption cross section. One can express the extinction cross sectionsin terms of extinction matrix elements

Cext =1

I inc

(K11(n

inc)I inc + K12(ninc)Qinc+

K13(ninc)U inc + K14(n

inc)V inc), (1.29)

and the scattering cross section in terms of phase matrix elements

Csca =1

I inc

∫4π

dr(Z11(r, ninc)I inc + Z12(r, ninc)Qinc+

Z13(r, ninc)U inc + Z14(r, ninc)V inc). (1.30)

The absorption cross section is the difference between extinction andscattering cross section:

Cabs = Cext − Csca. (1.31)

The single scattering albedo ω0, which is a commonly used quantityin radiative transfer theory, is defined as the ratio of the scatteringand the extinction cross section:

ω0 =Csca

Cext≤ 1. (1.32)

All cross sections are real-valued positive quantities and have the di-mension of area.

The phase function is generally defined as

p(r, ninc) =4π

CscaI inc

(Z11(r, ninc)I inc + Z12(r, ninc)Qinc+

Z13(r, ninc)U inc + Z14(r, ninc)V inc).(1.33)

The phase function is dimensionless and normalized:14π

∫4π

p(r, ninc) dr = 1. (1.34)

1.4 Particle Ensembles 27

1.4 Scattering, absorption and emissionby ensembles of independent particles

The formalism described in the previous chapter applies only for radia-tion scattered by a single body or a fixed cluster consisting of a limitednumber of components. In reality, one normally finds situations, whereradiation is scattered by a very large group of particles forming a con-stantly varying spatial configuration. Clouds of ice crystals or waterdroplets are a good example for such a situation. A particle collec-tion can be treated at each given moment as a fixed cluster, but as ameasurement takes a finite amount of time, one measures a statisticalaverage over a large number of different cluster realizations.

Solving the Maxwell equations for a whole cluster, like a collectionof particles in a cloud, is computationally too expensive. Fortunately,particles forming a random group can often be considered as inde-pendent scatterers. This approximation is valid under the followingassumptions:1. Each particle is in the far-field zone of all other particles.2. Scattering by the individual particles is incoherent.

As a consequence of assumption 2, the Stokes parameters of the par-tial waves can be added without regard to the phase. If the particlenumber density is sufficiently small, the single scattering approxima-tion can be applied. The scattered field in this approach is obtainedby summing up the fields generated by the individual particles in re-sponse to the external field in isolation from all other particles. Ifthe particle positions are random, one can show, that the phase ma-trix, the extinction matrix and the absorption vector are obtained bysumming up the respective characteristics of all constituent particles.

1.4.1 Single scattering approximation

We consider a volume element containing N particles. We assumethat N is sufficiently small, so that the mean distance between theparticles is much larger than the incident wavelength and the averageparticle size. Furthermore we assume that the contribution of the totalscattered signal of radiation scattered more than once is negligibly

28 1 Theoretical background

small. This is equivalent to the requirement

N 〈Csca〉l2

1, (1.35)

where 〈Csca〉 is the average scattering cross section per particle andl is the linear dimension of the volume element. The electric fieldscattered by the volume element can be written as the vector sum ofthe partial scattered fields scattered by the individual particles:

Esca(r) =N∑

n=1

Ensca(r). (1.36)

As we assume single scattering the partial scattered fields are givenaccording to Equation (1.17):(

[Escan (r)]θ

[Escan (r)]φ

)=

eikr

rS(r, ninc)

(Einc

0θ

Einc0φ

), (1.37)

where S is the total amplitude scattering matrix given by:

S(r, ninc) =N∑

n=1

ei∆nSn(r, ninc). (1.38)

Sn(r, ninc) are the individual amplitude matrices and the phase ∆n

is given by

∆n = krOn · (ninc − r), (1.39)

where the vector rOn connects the origin of the volume element O withthe nth particle origin (see Figure 1.2). Since ∆n vanishes in forwarddirection and the individual extinction matrices can be written interms of the individual amplitude matrix elements, the total extinctionmatrix is given by

K =N∑

n=1

Kn = N 〈K〉 , (1.40)

where 〈K〉 is the average extinction matrix per particle. One canderive the analog equation for the phase matrix

Z =N∑

n=1

Zn = N 〈Z〉 , (1.41)

1.4 Particle Ensembles 29

O

1

2

rO1

rO2

point

scattering medium

observation

r

Figure 1.2: A volume element of a scattering medium conststing of a particleensemble. O is the origin of the volume element, rO1 connects the originwith particle 1 and rO2 with particle 2. The observation point is assumedto be in the far-field zone of the volume element.

where 〈Z〉 is the average phase matrix per particle. In almost allpractical situations, radiation scattered by a collection of independentparticles is incoherent, as a minimal displacement of a particle or aslight change in the scattering geometry changes the phase differencesentirely. It is important to note, that the ensemble averaged phasematrix and the ensemble averaged extinction matrix have in general16 independent elements. The relations between the matrix elements,which can be derived for single particles, do not hold for particleensembles.

30 1 Theoretical background

1.5 Phenomenological derivation of theradiative transfer equation

When the scattering medium contains a very large number of particlesthe single scattering approximation is no longer valid. In this case wehave to take into account that each particle scatters radiation thathas already been scattered by another particle. This means that theradiation leaving the medium has a significant multiple scattered com-ponent. The observation point is assumed to be in the far-field zoneof each particle, but it is not necessarily in the far-field zone of thescattering medium as a whole. A traditional method in this case is tosolve the radiative transfer equation. This approach still assumes, thatthe particles forming the scattering medium are randomly positionedand widely separated and that the extinction and the phase matricesof each volume element can be obtained by incoherently adding the re-spective characteristics of the constituent particles. In other words thescattering media is assumed to consist of a large number of discrete,sparsely and randomly distributed particles and is treated as contin-uous and locally homogeneous. Radiative transfer theory is originallya phenomenological approach based on considering the transport ofenergy through a medium filled with a large number of particles andensuring energy conservation. Mishchenko (2002) has demonstratedthat it can be derived from electromagnetic theory of multiple wavescattering in discrete random media under certain simplifying assump-tions.

In the phenomenological radiative transfer theory, the concept ofsingle scattering by individual particles is replaced by the assumptionof scattering by a small homogeneous volume element. It is further-more assumed that the result of scattering is not the transformationof a plane incident wave into a spherical scattered wave, but the trans-formation of the specific intensity vector, which includes the Stokesvectors from all waves contributing to the electromagnetic radiationfield.

1.5 Radiative transfer equation 31

The vector radiative transfer equation (VRTE) is

dI(n, ν)ds

=− 〈K(n, ν, T )〉 I(n, ν) + 〈a(n, ν, T )〉B(ν, T )

+∫

4π

dn′ 〈Z(n,n′, ν, T )〉 I(n′, ν),(1.42)

where I is the specific intensity vector, 〈K〉 is the ensemble-averagedextinction matrix, 〈a〉 is the ensemble-averaged absorption vector, B isthe Planck function and 〈Z〉 is the ensemble-averaged phase matrix.Furthermore ν is the frequency of the radiation, T is the tempera-ture, ds is a path-length-element of the propagation path and n thepropagation direction. Equation (1.42) is valid for monochromatic orquasi-monochromatic radiative transfer. We can use this equation forsimulating microwave radiative transfer through the atmosphere, asthe scattering events do not change the frequency of the radiation.

The four-component specific intensity vector I = (I,Q,U, V )T fullydescribes the radiation and it can directly be associated with the mea-surements carried out by a radiometer used for remote sensing. Forthe definition of the components of the specific intensity vector referto Section 1.2, where the Stokes components are described. Since thespecific intensity vector is a superposition of Stokes vectors, the po-larization state of the specific intensity vector can be analysed in thesame way as the polarization state of the Stokes vector.

The three terms on the right hand side of Equation (1.42) de-scribe physical processes in an atmosphere containing different par-ticle types and different trace gases. The first term represents theextinction of radiation traveling through the scattering medium. Itis determined by the ensemble averaged extinction coefficient matrix〈K〉. For microwave radiation in cloudy atmospheres, extinction iscaused by gaseous absorption, particle absorption and particle scat-tering. Therefore 〈K〉 can be written as a sum of two matrices, theparticle extinction matrix 〈Kp〉 and the gaseous extinction matrix〈Kg〉:

〈K(n, ν, T )〉 = 〈Kp(n, ν, T )〉+ 〈Kg(n, ν, T )〉 . (1.43)

The particle extinction matrix is the sum over the individual specific

32 1 Theoretical background

extinction matrices 〈Kpi (n, ν, T )〉 of the N different particles types

contained in the scattering medium weighted by their particle numberdensities np

i :

〈Kp(n, ν, T )〉 =N∑

i=1

npi 〈K

pi (n, ν, T )〉 . (1.44)

A particle distribution, which can include various particle sizes, shapesand orientations, can be represented by a single particle type, since itis possible to derive an ensemble averaged phase matrix 〈Zi〉, an en-semble averaged extinction matrix 〈Ki〉 and an ensemble averaged ab-sorption vector 〈ai〉. The gaseous extinction matrix is directly derivedfrom the scalar gas absorption. As there is no polarization due to gasabsorption at cloud altitudes, the off-diagonal elements of the gaseousextinction matrix are zero. At very high altitudes above approximately40 km there is polarization due to the Zeeman effect, mainly due tooxygen molecules. However, in the toposphere and stratosphere molec-ular scattering can be neglected in the microwave frequency range.Hence the coefficients on the diagonal correspond to the gas absorp-tion coefficient:⟨

Kgl,m(ν, T )

⟩=

〈αg(ν, T )〉 if l = m

0 if l 6= m.(1.45)

where T is the temperature of the atmosphere and 〈αg〉 is the totalscalar gas absorption coefficient, which is calculated from the individ-ual absorption coefficients of all M trace gases αg

i (P, ν, T ) and theirvolume mixing ratios ng

i as:

〈αg(ν, T )〉 =M∑i=1

ngi α

gi (ν, T ). (1.46)

The second term in Equation (1.42) is the thermal source term. Itdescribes thermal emission by gases and particles in the atmosphere.The ensemble averaged absorption vector 〈a〉 is

〈a(n, ν, T )〉 = 〈ap(n, ν, T )〉+ 〈ag(ν, T )〉 , (1.47)

where 〈ap〉 and 〈ag〉 are the particle absorption vector and the gas

1.5 Radiative transfer equation 33

absorption vector, respectively. The particle absorption vector is asum over the individual absorption vectors 〈ap

i 〉, again weighted withnp

i :

〈ap(n, ν, T )〉 =N∑

i=1

npi 〈a

pi (n, ν, T )〉 . (1.48)

The gas absorption vector is simply

〈ag(ν, T )〉 = (〈αp(ν, T )〉 , 0, 0, 0)T . (1.49)

The last term in Equation (1.42) is the scattering source term. It addsthe amount of radiation which is scattered from all directions n′ intothe propagation direction n. The ensemble averaged phase matrix 〈Z〉is the sum of the individual phase matrices 〈Zi〉 weighted with np

i :

〈Z(n,n′, ν, T )〉 =N∑

i=1

npi 〈Zi(n,n′, ν, T )〉. (1.50)

The scalar radiative transfer equation (SRTE)

dI

ds(n, ν) = −〈K11(n, ν, T )〉 I(n, ν) + 〈a1(n, ν, T )〉B(ν, T )

+∫

4π

dn′ 〈Z11(n,n′, ν, T )〉 I(n′, ν) (1.51)

can be used presuming that the radiation field is unpolarized. This ap-proximation is reasonable if the scattering medium consists of spheri-cal or completely randomly oriented particles, where 〈Kp〉 is diagonaland only the first element of 〈ap〉 is non-zero.

2 ARTS – the atmosphericradiative transfer system

This chapter introduces basic concepts and definitions of the ARTSmodel. It provides a brief summary of functions and methods usedfor the scattering simulations. Many of the functions, for examplefunctions for the calculation of propagation paths, could be sharedbetween the clear sky part of the model and the scattering part.

2.1 HistoryA lot of effort has been put in the development of dedicated forwardmodels for different sensors. All of these models have many parts incommon. While appropriate for operational data analysis, such spe-cialized models are not appropriate for scientific studies of new sensorconcepts, since they can not easily be adapted to new instruments.This was the reason for the development of more general forwardmodels like the program FORWARD (Eriksson and Buehler, 2001),which was mostly written by J. Langen in the time period 1991–1998at the University of Bremen, or the Skuld model mainly developedby Eriksson et al. (2002) during 1997–1998. Although these modelswere rather general and have been used successfully over the years,both suffered from being not easily modifiable and extendable. Thishas lead to the development of a model which emphasizes modularity,extendibility, and generality.

It was decided that the development work should be shared betweenthe Bremen and Chalmers universities, with Bremen being largelyresponsible for the overall program architecture and the absorptionpart, Chalmers being largely responsible for the radiative transfer part

35

36 2 ARTS

and the calculation of Jacobians. The project was put under a GNUgeneral public license (Stallman, 2002), in order to give the right legalframework for such a true collaboration.

The program, along with extensive documentation, is freely avail-able on the Internet, under http://www.sat.uni-bremen.de/arts/. Thestable 1-0-x branch of the program is described in Buehler et al.(2005a). Stable means that there will be only bug fixes, no additionsof new features.A great part of the work for this thesis is dedicatedto the development of the new branch, 1-1-x, which can handle scat-tering in the atmosphere. The development team has been joined byC. Davis from the University of Edinburgh, who has used the ARTSmodel as a platform for the implementation of a Monte Carlo scat-tering module additionally to the discrete ordinate scattering module,which is presented in this thesis.

2.2 Definition of the atmosphere2.2.1 Atmospheric dimensionality

The modeled atmosphere can be selected to have different dimension-alities:

3D This is the most general case, where the atmospheric fields varyin all three spatial coordinates. A spherical coordinate system isused where the dimensions are pressure (P ), latitude (α) and longi-tude (β). Choosing this option allows to simulate realistic radiationfields, including strongly horizontally inhomogeneous cloud cover-age.

1D A “1D” atmosphere is a spherically symmetric atmosphere, whichmeans that atmospheric fields and the ground extend in all threedimensions, but they do not have a variation in latitude and lon-gitude. Atmospheric fields for instance vary only as a function ofaltitude. The surface of the earth corresponds to a sphere. The 1Dgeometry is a crude approximation for scattering simulations, asa spherically symmetric cloud corresponds to a globally complete

2.2 Definition of the atmosphere 37

cloud coverage. This extreme case can be used to study the effectof scattering by largely extended thin cirrus clouds.

2D A 2D atmosphere extends inside a plane. A polar coordinatesystem, consisting of a radial and an angular coordinate, is used.The 2D case is most likely used for satellite measurements wherethe atmosphere is observed inside the orbit plane. Scattering cal-culations can not be performed in 2D geometry as there is no caseinvolving clouds that give rise to a radiation field that fits into a2D framework.

2.2.2 The cloud box

In order to save computational time, scattering calculations are lim-ited to the part of the atmosphere containing clouds and other scat-tering objects. The atmospheric region in which scattering shall beconsidered is denoted as the cloud box. The cloud box is defined tobe rectangular in the used coordinate system, with limits exactly atpoints of the involved grids. This means, for example, that the verticallimits of the cloud box are two pressure surfaces.

When defining the cloud box limits, one must avoid that radiationemerging from the cloud box reenters the cloud box at another point,because the scattering calculation takes as boundary condition theincoming clear sky field. If there is a large amount of ground reflectionand if the optical depth under the cloud box is small, there shouldnot be a gap between the ground and the cloud box. In this case theground is included as a scattering object. However, for many cases itcan be accepted to have a gap between the ground and the cloud box,with the gain that the cloud box can be made smaller. Such a case iswhen the ground is treated to act as blackbody, the ground is thennot reflecting any radiation. Reflections from the ground can also beneglected if the zenith optical thickness of the atmosphere betweenthe ground and cloud box is sufficiently high.

Figure 2.1 shows schematically 1D and 3D model atmospheres in-cluding a cloud box.

38 2 ARTS

Atmospheric gridsGroundGeoidCloud boxAtmospheric field

1D

3D

Figure 2.1: Top panel: Schematic of a 1D atmosphere. The atmosphere ishere spherically symmetric. This means that the radius of the geoid, theground and all atmospheric profiles are constant around the globe. Thegrey area indicates the cloud box. Bottom panel: Cross section of a 3Datmosphere. Atmospheric fields are defined on all grid points. The cloud boxhas a finite horizontal extent and the surface is not spherically symmetric.

2.3 Radiative transfer calculationsThe radiative transfer (RT) calculations are divided into two separateparts, a clear sky part and the scattering calculations inside the cloudbox. These parts have been implemented as two main modules with awell defined interface. The task of the scattering part is to determinethe outgoing intensity field of the cloud box. The scattering calcula-tions can be performed in any way as long as the outgoing field isprovided. The outgoing field is then used as the radiative backgroundfor observation directions giving a propagation path that intersectswith the cloud box.

The aim, when designing the clear sky and scattering modules, wasthat as many components as possible should be common. This is ad-vantageous for many reasons, e.g., it decreases the amount of code

2.3 Radiative transfer calculations 39

to maintain, it facilitates detection of bugs, and enhances the con-sistency between the modules. An example of a common componentis the calculation of propagation paths, where the same function forpropagation path calculations is used in both modules.

The clear sky radiative transfer calculation is performed for a fullmeasurement sequence. This means that the function calculates spec-tra for all positions of the sensor, and all pencil beam directions neededfor the weighting with the antenna pattern. The inclusion of sensorcharacteristics etc. are not discussed here, details are given in Erikssonet al. (2004).

The clear sky part is vectorized in frequency (all monochromaticfrequencies are handled in parallel), being in contrast to the scatteringpart. The number of Stokes components to consider can be set to anyvalue from 1 to 4 (this is also valid for the scattering part). Polarizedcalculations (number of Stokes components > 1) can be performedindependently from the cloud box being activated or not.

The RT calculations for a single pencil beam direction can be sep-arated into three sub-tasks:

– Calculation of the propagation path.– Determining the radiative background.– Solving the radiative transfer equation.

2.3.1 Propagation paths

Any combination of sensor position and line-of-sight that makes sensewith respect to the model atmosphere is allowed. The main restrictionis that propagation paths are only allowed to enter or exit the modelatmosphere at the top. This means that the propagation path can notexit the model atmosphere at a latitude end face for a 2D and 3Dcase.

Propagation paths are calculated backwards from the sensor to thepractical starting point. If the sensor is placed outside the model atmo-sphere, geometrical calculations are used to find the exit point at thetop of the atmosphere. Inside the model atmosphere, the path is cal-culated in steps, from one crossing of a grid cell boundary to the next

40 2 ARTS

Figure 2.2: Propagation path examples for a 1D atmosphere. The dottedlines are the atmospheric grid, the dashed-dotted line is the geoid, the thicksolid line is the ground and the cylindrical segment drawn with a thin solidline is the cloud box.

crossing. The functions to calculate such propagation path steps canbe used both in the clear sky and scattering parts. Propagation pathsare followed backwards until the top of the atmosphere, the groundor the cloud box (if activated) is reached. The propagation paths aredescribed by a number of points along the path. Points are always in-cluded for crossings with the grids, tangent points (if such exist) andthe position of the sensor, if placed inside the atmosphere. Dependingon the function selected for the path step calculations, other pointscan be included, for instance to fulfill a criterion on the maximumlength along the path between the points. Examples of propagationpaths are given in Figure 2.2.

2.3.2 Radiative background

The radiative intensities at the starting point of the propagation pathare denoted as the radiative background. Four possible radiative back-grounds exist:

1. Cosmic background when the propagation path starts at the top ofthe atmosphere.

2. The up-welling radiation from the ground when the propagation

2.3 Radiative transfer calculations 41

path intersects with the ground. Emission and scattering proper-ties of the ground need to be defined to determine the radiativebackground.

3. If the propagation path hits the surface of the cloud box, the ra-diative background is obtained by interpolating the radiation fieldleaving the cloud box. This interpolation considers the propagationdirection of the path at the crossing point.

4. The internal intensity field of the cloud box is the radiative back-ground for cases when the sensor is placed inside the cloud box.

2.3.3 Clear sky radiative transfer

The intensity matrix (holding all frequencies and Stokes components)is set to equal the radiative background. The calculations are thenperformed by solving the radiative transfer problem from one point ofthe propagation path to next, until the end point is reached.

The clear sky vector radiative transfer equation follows from thegeneral VRTE Equation (1.42) by omitting the scattering integraland particle contributions to the extinction matrix and the absorptionvector:

dI

ds(n, ν, T ) = −〈Kg(n, ν, T )〉I(n, ν, T ) + 〈ag(n, ν, T )〉B(ν, T ).

(2.1)

This equation can be solved analytically for constant coefficients. Theextinction matrix 〈Kg(n, ν, T )〉 and the absorption vector 〈ag(n, ν, T )〉are averaged for one propagation path step. The averaging procedurewill be described more detailed in Section 4.1. The solution is foundusing a matrix exponential approach (see Appendix B.1):

Ii = e−〈Kg〉s · Ii−1 + (I− e−〈K

g〉s)〈Kg〉−1

(〈ag〉 B), (2.2)

where 〈Kg〉 and 〈ag〉 are the averaged quantities and i denotes a pointin the propagation path.

42 2 ARTS

2.4 ScatteringAs mentioned above the task of the scattering module is to deter-mine the outgoing radiation field on the boundary of the cloud box.This requires numerical methods to solve the VRTE (1.42) inside thecloud box as there is no analytical solution to the VRTE without anyapproximations. Two different approaches are implemented in ARTS:A backward Monte Carlo scheme which is briefly described in Sec-tion 5.3 and the discrete ordinate iterative approach, which has beendeveloped by the author of this thesis and will be described in detailin Chapter 4. Several studies in which the DOIT method has beenapplied will be presented in the chapters 6 to 8.

2.5 Gas absorptionCalculating gas absorption in a line-by-line way is expensive, as some-times contributions from thousands or ten thousands of spectral lineshave to be taken into account. This needs to be done over and overagain for each point in the atmosphere. The gas absorption coefficientdoes not depend directly on the position, but on the atmospheric statevariables: pressure, temperature and trace gas concentrations. The ba-sic idea in ARTS is to pre-calculate absorption for discrete combina-tions of these variables, store the values in a lookup table, and theninterpolate them for the actual atmospheric state. The gas absorptioncoefficients are taken from spectral line catalog, for example from theHITRAN catalogs (Rothman et al., 1998).

2.6 Definition of clouds and atmosphericfields

In the Earth’s atmosphere we find liquid water clouds consisting ofapproximately spherical water droplets and cirrus clouds consisting ofice particles of diverse shapes and sizes. We also find different kinds ofaerosols. In order to take into account this variety, the model allows to

2.7 Unit conventions 43

define several particle types. A particle type is either a specified parti-cle or a specified particle distribution, for example a particle ensemblefollowing a gamma size distribution. The particles can be completelyrandomly oriented, azimuthally randomly oriented or arbitrarily ori-ented. For each particle type being a part of the modeled cloud field,a data file containing the single scattering properties (〈Ki〉, 〈ai〉, and〈Zi〉), and the appropriate particle number density field is required.The particle number density fields are stored in data files, which in-clude the field stored in a three-dimensional tensor and also the ap-propriate atmospheric grids (pressure, latitude and longitude grid).For each grid point in the cloud box the single scattering propertiesare averaged using the particle number density fields. In the scatteringdatabase the single scattering properties are not always stored in thesame coordinate system. For instance for randomly oriented particlesit makes sense to store the single scattering properties in the so-calledscattering frame in order to reduce memory requirements (refer toSection 3.4 for more details).

The atmospheric fields, which are temperature, altitude, and vol-ume mixing ratio fields, are stored in the same format as the particlenumber density fields.

2.7 Unit conventionsInternally the ARTS model uses SI1 units for all quantities. However,SI units can sometimes be inconvenient, for example to represent theradiation field in the atmosphere. Therefore it is possible in ARTS toconvert radiances from the SI unit [W s m−2 sr−1] into a brightnesstemperature [K] unit. There are two brightness temperature (BT)definitions which can be applied:1. Planck BT:

The simplest case of remote sensing of temperature occurs whenatmospheric extinction can be neglected. Then a satellite would just’see’ the thermal emission of the Earth’s surface. One obtains the

1 SI units – Système International d’Unités

44 2 ARTS

temperature of the surface from the measured radiance by invertingthe Planck function Equation (1.26):

TPlanck(Ib) =hν

kB ln(

2hν3

c2Ib+ 1) . (2.3)

More generally we can define a brightness temperature in terms ofthe radiance Ib using Equation (2.3), even in presence of extinction.

2. Rayleigh Jeans BT:The definition of Planck BT can not be used for all Stokes compo-nents, because it is only defined for positive values of Ib. Since theStokes components Q, U and V can be negative, we would like tohave a conversion, which can also be applied for negative values.Another problem is the non-linearity of the Planck BT definition.The Stokes component Q is the difference between the vertically andthe horizontally polarized parts of the radiation. Using the PlanckBT definition, the value of Q would depend on whether we firsttransform Iv and Ih and take the difference of the obtained PlanckBT, or we transform Q directly to Planck BT. The Rayleigh Jeansdefinition of BT is linear. The proportionality factor is derived fromthe Rayleigh Jeans approximation: At small frequencies (hν kbT )the Planck function is approximately

Ib(T, ω) =2kBν2

c2T. (2.4)

Inverting Equation (2.4) yields the definition of Rayleigh Jeans BT:

TRJ(Ib) =c2

2kBν2Ib. (2.5)

Since Rayleigh Jeans BT’s can be used for all Stokes components,this BT definition is the only one used in this work. Note: RayleighJeans BT are not equal to Planck BT, they are two different unitsto represent radiances.

3 Description of cloudsas scattering media

This chapter deals with the representation of clouds as scattering me-dia in radiative transfer models. An overview of cloud microphysicsprovides realistic ranges of particle sizes, shapes and ice mass contentsof cirrus clouds. Different methods to calculate scattering propertiesfor small particles are introduced. It is shown that the Rayleigh andMie approximations are not sufficient for modeling radiative transferthrough cirrus clouds. The cloud particles are not sufficiently smallto be treated as Rayleigh scatterers. They are usually aspherical andoften horizontally aligned, which makes it impossible to use the Mietheory, which is valid only for spherical particles. Although it canonly handle rotationally symmetric particles, the T-matrix methodwas chosen to be used in ARTS, since it yields a rather good approx-imation for most realistic particles and it is widely used and tested.This chapter also introduces particle size distributions which are usedfor simulations in later chapters.

3.1 Microphysics of cloudsThe earth’s atmosphere consists of various particles: aerosols, waterdroplets, ice crystals, raindrops, snowflakes, graupel and hailstones.Cloud particles are the most important scatterers in the microwaveregion.

Clouds, which are composed of water droplets or ice crystals, areconventionally classified in terms of their position and appearance inthe atmosphere. At mid-latitudes, clouds with base heights of about6 km are defined as high clouds or cirrus clouds. The group of low

45

46 3 Clouds as scattering media

clouds below about 2 km include stratus and cumulus. Middle clouds,between the high and the low clouds, include altocumulus and al-tostratus. The dispersion of particle sizes and their phase (liquid orice) determines the microphysical state of a cloud. According to Liou(2002) low clouds and some middle clouds are generally composedof spherical water droplets with sizes ranging from 1 µm to 20 µm.The typical size for a water droplet is 5 µm. Middle clouds with tem-peratures warmer than about −20C can contain super-cooled waterdroplets that coexist with ice particles. The small water droplets arespherical due to the surface tension. Larger raindrops deviate fromthe spherical shape while they are falling down. Cirrus clouds andsome of the top and middle clouds contain ice crystals. The ice crys-tal shapes are irregular and depend on temperature, relative humidityand on the dynamics in the clouds, i.e., whether they undergo collisionand coalescence processes. For humidities close to water saturation,the particles have prismatic skeleton shapes that occur in hollow andcluster crystals. These particles are referred to as bullet rosettes andthey occur for example in cirrocumulus clouds. In cirrostratus clouds,where the relative humidity is close to ice saturation, ice crystals arepredominantly individual and have shapes like columns, prisms, andplates. Between water and ice saturation, the ice crystals grow in theform of prisms.

Figure 3.1 shows a spectrum of ice crystal sizes and shapes as afunction of height, relative humidity, and temperature in a typicalmidlatitude cirrus. Since ice crystal shape and size vary greatly withtime and space, it is difficult to find representative values for remotesensing applications. Figure 3.2 shows five measured size distributions.The data is taken from Heymsfield and Platt (1984) and the figure isadapted from Liou (2002). The mean effective ice crystal size rangesfrom 10 µm to 124 µm.

In cirrus clouds ice particles are generally not randomly oriented.Laboratory experiments have shown, that cylinders (aspect ratio1 <

1) tend to fall with their long axes horizontally oriented. Based onobservations, columnar and plate crystals (aspect ratio > 1) tend to

1 The aspect ratio of a particle is its diameter divided by its length.

3.1 Microphysics of clouds 47

Figure 3.1: Ice crystal shape and size as a function of height and relativehumidity captured by a replicator ballon sounding system in Marshall, Col-orado on November 10, 1994. The relative humidity was measured by acryogenic hygrometer (dashed line) and Vaisala RS80 instruments (solidline and dots). Figure adapted from Liou (2002).

fall with their major axes horizontally oriented. The orientation of iceparticles in cirrus clouds has been observed by lidar measurementsbased on the depolarization technique in the backscattering direction(Liou, 2002). The measurements have shown, that specific orientationsoccur when the particles have relatively large sizes and well definedshapes, like cylinders or plates. If the ice crystals are irregular, suchas aggregates, there is no preferred orientation. Moreover, smaller iceparticles in cirrus clouds tend to be not randomly oriented in three-dimensional space if there is a substantial turbulence in the cloud. Ithas also been observed that ice particle orientation and alignment arestrongly modulated by the electric field in clouds.

48 3 Clouds as scattering media

Figure 3.2: Ice crystal size distributions for midlatitude cirrus clouds cover-ing a range of mean effective ice crystal sizes from 10 µm (Contrail), 24 µm(Cold), 42 µm (Cirrostratus), 75 µm (Thick), to 124 µm (Uncinus). The datawas taken from Heymsfield and Platt (1984) and from Liou et al. (1998).Figure adapted from Liou (2002).

3.2 Coordinate systems: The laboratoryframe and the scattering frame

For radiative transfer calculations we need a coordinate system to de-scribe the direction of propagation. For this purpose we use the labo-ratory frame, which has been introduced in Chapter 1, Figure 1.1. Thez-axis corresponds to the local zenith direction and the x-axis pointstowards the north-pole. The propagation direction is described by thelocal zenith angle θ and the local azimuth angle φ. This coordinate

3.2 Coordinate systems 49

incn

nsca

x

z

θ

Figure 3.3: Illustration of the scattering frame. The z-axis coincides withthe incident direction ninc. The scattering angle Θ is the angle between ninc

and nsca.

system is the most appropriate frame to describe the propagation di-rection and the polarization state of the radiation. However, in orderto describe scattering of radiation by a particle or a particle ensemble,it makes sense to define another coordinate system taking into con-sideration the symmetries of the particle or the scattering medium, asone gets much simpler expressions for the single scattering properties.For macroscopically isotropic and mirror-symmetric scattering mediait is convenient to use the scattering frame, in which the incidencedirection is parallel to the z-axis and the x-axis coincides with thescattering plane, that is, the plane through the unit vectors ninc andnsca. The scattering frame is illustrated in Figure 3.3. For symmetryreasons the single scattering properties defined with respect to thescattering frame can only depend on the scattering angle Θ,

Θ = arccos(ninc · nsca), (3.1)

between the incident and the scattering direction.

50 3 Clouds as scattering media

3.3 Methods to calculate scattering bysmall particles

This section introduces basic concepts, which are relevant for atmo-spheric scattering. It will be shown that neither the Rayleigh scatter-ing approximation nor the Lorentz-Mie theory for scattering of radia-tion by small particles is generally appropriate to describe scattering ofmicrowave radiation by cirrus cloud particles. The T-matrix method,one of the more elaborate methods, leads to a better understanding,especially of polarization due to cloud scattering.

An important quantity in scattering theory for small particles is thesize parameter

x =2πr

λ=

2πrν

c, (3.2)

where r is the particle radius, λ is the wavelength and ν is the fre-quency of the incident radiation. The size of non-spherical particlesis here defined by their equal volume sphere radius. The volume ofa cylindrical particle for instance, which has an equal volume sphereradius of 75 µm, is identical to the volume of a sphere with a radius of75 µm. Rayleigh scattering, the most simple theory, is valid for x 1,i.e., if the particle size is much smaller than the wavelength of theincident radiation. If the wavelength is comparable to the particle size(x ≈ 1), one can apply the Lorentz-Mie theory for spherical particlesor the T-matrix method for spherical and non-spherical particles. Forsize parameters much greater than one, the geometrical optics approx-imation can be applied. Figure 3.4 shows size parameters as a functionof frequency for different particle sizes of cloud ice particles. Only forvery small cloud particles (10 µm) or frequencies below 100 GHz x issmall enough such that Rayleigh scattering applies. The geometricaloptics approximation can not be applied for any of these particles.

3.3.1 Rayleigh scattering

Rayleigh scattering is mostly applied for molecular scattering in thevisible wavelength region or for the scattering of very low-frequency

3.3 Computation methods 51

10 100 10000.001

0.01

0.1

1

10

Frequency [ GHz ]

x [ ]

10 µm

42 µm

75 µm

124 µm

Figure 3.4: Size parameter x as a function of frequency for different particlesizes typical for cloud particles.

microwave radiation by hydrometeors. As air molecules are of severalorders smaller than microwave wavelengths, molecular scattering canbe neglected in the microwave wavelength region. Nevertheless, someaerosol or cloud particles can be sufficiently small to be treated asRayleigh scatterers. The classical electro-dynamical solution yields therelation between incident and scattered intensity, which are denotedby I0 and I respectively,

I =I0

r2p

α2

(2π

λ

)4 1 + cos2 Θ2

, (3.3)

where α is the polarizability of the small particle, rp is the distancefrom the particle and Θ is the scattering angle. This is the formuladerived by Rayleigh (1871). The formula shows that the intensity ofsunlight scattered by a molecule is proportional to the incident inten-sity and inversely proportional to the square of the distance betweenthe molecule and the observation point. It is also inversely propor-tional to λ4 which explains the blue color of the sky. Since blue lightis scattered more than red light, the sky, when viewed away from thesun disk appears blue. The Rayleigh phase function p(Θ) is given by

p(Θ) =34(1 + cos2 Θ). (3.4)

The phase function is symmetric about the minimum at a scatteringangle of 90. An equal amount of radiation is scattered into the for-

52 3 Clouds as scattering media

ward direction (Θ = 0) and into the backward direction (Θ = 180).The degree of linear polarization plin (Equation (1.15)) can also be de-rived for particles which are very small compared to the wavelength:

plin =sin2 Θ

cos2 Θ + 1. (3.5)

In the forward and backward directions the scattered radiation re-mains completely unpolarized, whereas at a scattering angle of 90,the scattered radiation becomes completely polarized. In other direc-tions the radiation is partially polarized.

3.3.2 Lorentz-Mie theory for scatteringby spherical particles

Starting from the Maxwell equations, one can derive analytically thesingle scattering properties of spherical particles in the far field ap-proximation. For a detailed derivation refer, for example, to Liou et al.(1998) or Bohren and Huffman (1998). The MATLAB code developedby Mätzler (2002) was used to calculate phase functions for differentparticle sizes existing in clouds, corresponding to the effective radii ofthe size distributions shown in Figure 3.2. The particles were assumedto consist of ice. The refractive index of ice is calculated according toMätzler (1998). The phase functions are computed for 89 GHz and318 GHz and presented in Figure 3.5. For 89 GHz all phase functionsare very close to the Rayleigh phase function, which is the expectedresult, as the size parameter for this frequency is always smaller thanapproximately 0.01. However, for 318 GHz the phase functions are verydifferent from the Rayleigh phase function. Only for a particle radiusof 10 µm it is still possible to use the Rayleigh approximation. Thelarger particles scatter more radiation into the forward and less intothe backward direction. The minimum of the Rayleigh phase functionis at a scattering angle of 90. It shifts towards larger scattering anglesas the size parameter increases.

Figure 3.6 shows the extinction efficiency Qext, the scattering effi-ciency Qsca and the absorption efficiency Qabs of spherical ice particles

3.3 Computation methods 53

0 20 40 60 80 100 120 140 160 1800.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

1.6

Theta [ °]

P [

]89 GHz

0 20 40 60 80 100 120 140 160 180

0.8

1

1.2

1.4

1.6

1.8

2

2.2

Theta [ ° ]

P [

]

318 GHz

Rayleigh

10 µm

42 µm

75 µm

124 µm

Figure 3.5: Grey lines: Rayleigh phase function. Black lines: Mie phase func-tions for different particle sizes. The left plot is for 89 GHz and the rightplot for 318 GHz.

in the whole microwave-wavelength region. The efficiencies correspondto the area normalized optical cross sections Cext, Csca and Cabs, whichhave been defined in Section 1.3.5, thus

Qext =Cext

πr2, Qsca =

Csca

πr2, Qabs =

Cabs

πr2, (3.6)

where r is the radius of the particle. For frequencies below 1000 GHz,the absorption efficiency is at least one order of magnitude smallerthan the scattering efficiency. Therefore the extinction efficiency andthe scattering efficiency are almost identical in this frequency range.The major maxima and minima of the scattering efficiency are calledthe interference structure and the irregular structure is called theripple structure. The origin of the term interference structure lies inthe interpretation of extinction as the interference between the inci-dent and the forward-scattered light. The scattering efficiency Qsca

increases rapidly until the size parameter reaches approximately twofor non-absorbing ice-particles. This means that for larger particlesthe maximum shifts towards lower frequencies. In later chapters, sim-ulations for satellite limb measurements at 318 GHz will be shown.Figure 3.6 shows, that the scattering signal at this frequency shoulddepend very much on the particle size. Another result of the Lorentz-Mie theory is that for a non-absorbing medium, for which the imagi-nary part of the refractive index equals zero, the scattering efficiency

54 3 Clouds as scattering media

approaches an asymptotic value of two for large size parameters. Thisimplies that the particle removes exactly twice the amount of en-ergy that it can intercept. It includes the diffracted component, whichpasses by the particle, and additionally the radiation scattered byreflection and refraction inside the particle.

Figure 3.7 shows the extinction efficiency Qext, the scattering ef-ficiency Qsca and the absorption efficiency Qabs of spherical waterdroplets of typical sizes. In contrast to ice particles, liquid waterdroplets mainly absorb microwave radiation. Below 1000 GHz scat-tering is negligibly small.

For polarized radiative transfer calculations, phase function, extinc-tion coefficient, scattering coefficient and absorption coefficient are notsufficient. The VRTE Equation (1.42) shows that we need the phasematrix 〈Z〉, the extinction matrix 〈K〉 and the absorption coefficientvector 〈a〉. The phase matrix represented in the scattering frame iscommonly called scattering matrix F . It follows from the Mie theorythat the scattering matrix has only four independent matrix elements;it reduces to the simple form:

F (Θ) =

F11(Θ) F12(Θ) 0 0F12(Θ) F11(Θ) 0 0

0 0 F33(Θ) F34(Θ)0 0 −F34(Θ) F33(Θ)

. (3.7)

For spherical particles F depends only on the scattering angle Θ, whichis obvious for symmetry reasons.

3.3.3 T-matrix method

The deficit of the Lorentz-Mie theory is, that it provides scatteringproperties only for spherical particles. As shown in Section 3.1, iceparticles are usually not spherical. The T-matrix method, which wasinitially introduced by Waterman (1965), can be applied for the com-putation of electromagnetic scattering by single, homogeneous, arbi-trarily shaped particles. This original method is also known as theextended boundary condition method (EBCM). At present, the T-matrix approach is one of the most powerful and widely used tool for

3.3 Computation methods 55

0 500 1000 1500 2000 2500 30000

1

2

3

4

5Q

ext

0 500 1000 1500 2000 2500 30000

1

2

3

4

5

Qsc

a

500 1000 1500 2000 2500

0.2

0.4

0.6

0.8

1

1.2

frequency [ GHz ]

Qab

s

42 µm75 µm124 µm248 µm

Figure 3.6: Extinction efficiency Qext, scattering efficiency Qsca and ab-sorption efficiency Qabs in the microwave-wavelength region for sphericalice particles with radii of 42 µm, 75 µm, 124 µm and 248 µm.

56 3 Clouds as scattering media

0 500 1000 1500 2000 2500 30000

0.2

0.4

0.6

0.8

Qex

t

0 500 1000 1500 2000 2500 30000

0.1

0.2

0.3

0.4

Qsc

a

1 µm5 µm10 µm20 µm

500 1000 1500 2000 2500

0.1

0.2

0.3

0.4

frequency [ GHz ]

Qab

s

Figure 3.7: Extinction efficiency Qext, scattering efficiency Qsca and absorp-tion efficiency Qabs in the microwave-wavelength region for spherical liquidwater droplets with radii of 1 µm, 5 µm, 10 µm and 20 µm.

3.3 Computation methods 57

the computation of scattering by aspherical single and compoundedparticles. A detailed theoretical explanation can be found in the bookby Mishchenko et al. (2002), who has developed several public domainprograms which are available at http://www.giss.nasa.gov/∼crmim.The T-matrix method needs the refractive index of ice as an input.Here the data from Warren (1984) was used to obtain the refractiveindex for the required frequencies and temperatures.

T-matrix program for randomly oriented, homogeneous androtationally symmetric particles

The program for randomly oriented rotationally symmetric particlescan be used for mono-disperse particles or for several analytical sizedistributions, e.g., the gamma distribution.

The total scattering matrix for the particle distribution is

F (Θ) =

F11(Θ) F12(Θ) 0 0F12(Θ) F22(Θ) 0 0

0 0 F33(Θ) F34(Θ)0 0 −F34(Θ) F44(Θ)

= N 〈F (Θ)〉 ,

(3.8)

where N is the number of particles in a unit volume and 〈F (Θ)〉 isthe ensemble-averaged scattering matrix per particle. In contrast tothe scattering matrix for spherical particles, the matrix elements F11

and F22 as well as F33 and F44 are different, therefore we have nowsix instead of four independent elements. Like for spherical particles,the scattering matrix depends only on the scattering angle Θ.

The T-matrix program allows the computation of the scatteringproperties for non-spherical, rotationally symmetric particles. Theshape of an aspherical particle is defined by its aspect ratio, whichis the diameter of the particle divided by its length. Thus, a particlewith an aspect ratio larger than one is a oblate particle and a particlewith an aspect ratio smaller than one is a prolate particle.

Single scattering properties for mono-disperse spheroidal particleswith an equal volume sphere radius of 75 µm are calculated. The re-sults of the calculations for 318 GHz are shown in Figure 3.8. All the

58 3 Clouds as scattering media

six phase matrix elements are presented. They are normalized by mul-tiplication with 4π

〈Csca〉 , where 〈Csca〉 is the ensemble averaged scatter-ing cross-section. The grey lines correspond to the Mie calculation andthe black lines correspond to the T-matrix calculations. The solid linescorrespond to the normalized scattering matrix for spherical particles,which are identical to randomly oriented spheroids with an aspect ra-tio of 1.0. The dashed lines correspond to prolate spheroids with aspectratios of 0.3 (thick) and 0.5 (thin). The dotted lines are the results foroblate spheroids with aspect ratios of 5.0 (thick) and 2.0 (thin). TheMie result corresponds to the T-matrix result for spherical particles,which means that the two methods are consistent. The matrix elementF11 which corresponds to the phase function, shows that randomlyoriented aspherical particles scatter slightly more radiation into theforward and slightly less radiation into the backward direction com-pared to spherical particles. The difference increases with increasingdeformation. The matrix element F21 = F12 mainly determines thepolarization state of the scattered radiation, since, for an unpolarizedincident beam, the second Stokes component Q of the scattered beamcorresponds to the product of F21 and the incoming intensity I0. Theplot shows, that maximal polarization occurs at a scattering angle ofabout 90ř and that Q is negative. The matrix element F22, whichequals F11 for spherical particles, is for aspherical particles smallerthan F11, especially in the backward direction. The absolute valueof the element F34 = F43 is very small. F33 and F44 deviate onlyslightly from the result for spherical particles. Overall, the calculationsfor randomly oriented aspherical particles show, that at 318 GHz fora particle size of 75 µm the phase matrix is very similar to that forspherical particles. Therefore the dependence of simulated radianceson the particle shape at this frequency is expected to be much lessthan the dependence on the particle size.

T-matrix program for a particle in arbitrary orientation

All previous results were obtained for randomly oriented particles.In reality the ice crystals in cirrus clouds tend to be horizontallyaligned. Another T-matrix program by Mishchenko (2000) is appli-

3.3 Computation methods 59

0 50 100 1500.5

1

1.5

2

Θ [ ° ]

F11

[ ]

Mie aspect ratios: 0.3 0.5 1.0 2.0 5.0

0 50 100 150−0.8

−0.6

−0.4

−0.2

0

Θ [ ° ]

F12

[ ]

0 50 100 1500.5

1

1.5

2

Θ [ ° ]

F22

[ ]

0 50 100 150−2

−1

0

1

2

Θ [ ° ]

F33

[ ]

0 50 100 150−1.5

−1

−0.5

0x 10−3

Θ [ ° ]

F34

[ ]

0 50 100 150−2

−1

0

1

2

Θ [ ° ]

F44

[ ]

Figure 3.8: Grey line: Phase matrix elements calculated using Lorentz-Mie theory. Black lines: Orientation averaged phase matrix elements forspheroidal particles with different aspect ratios calculated using the T-matrix program for randomly oriented particles. The dashed lines corre-spond to prolate particles with aspect ratios of 0.3 (thick) and 0.5 (thin).The dotted lines correspond to oblate particles with aspect ratios of 2.0(thin) and 5.0 (thick). The solid line corresponds to spherical particles. Theequal volume sphere radius is 75 µm and the frequency is 318 GHz.

60 3 Clouds as scattering media

cable for aspherical rotationally symmetric particles in arbitrary ori-entation. This program has been used to calculate the phase matrixfor cylindrical particles with different aspect ratios, which are hori-zontally oriented. The equal volume sphere radius is again 75 µm andthe frequency of the calculation is 318 GHz. To be able to comparethe results with Mie calculations, we have calculated the phase ma-trix for θinc = 0, φinc = 0, and φsca = 0, so that it corresponds to thescattering matrix F ,

F (θsca) = Z(θsca, φsca = 0, θinc = 0, φinc = 0). (3.9)

The z-axes of the laboratory frame, in which the phase matrix iscalculated, is chosen to be parallel to the incident direction (θinc =0, φinc = 0). Thus θsca corresponds to the scattering angle Θ. The par-ticle frame is defined in such a way that the z′-axes is parallel to thesymmetry axes of the particle. Figure 3.9 shows that for a horizontallyoriented plate the laboratory frame corresponds to the particle frame.For the horizontally oriented cylinder the z′-axes of the particle frameis rotated about the x-axes by 90 w.r.t. the z-axes. The orientation ofthe particles can be specified by appropriate Euler angles of rotation.The results of the calculations are presented in Figure 3.10. A cylin-der with an aspect ratio of 1.0 gives almost identical results to theMie calculation, when the symmetry axes corresponds to the z-axes.The results for the plates (dotted lines), which are oriented horizon-tally, are still similar. But the results for oriented cylinders (dashedlines) are very different. This is reasonable considering the symmetry.Like a sphere the plate has a circular geometrical cross-section whenit is seen from the top. The cylinder, which is oriented horizontally,has a rectangular geometrical cross-section when seen from the top.The linear polarization is the difference between the vertical and thehorizontal component of the intensity. As the horizontal and the ver-tical component of the electric field vector are perpendicular to thedirection of propagation, the plate must have the same influence onboth components for forward and backward scattering, because it issymmetric about the propagation direction of the radiation. ThereforeF21, which is related to linear polarization, must be zero for Θ = 0

and for Θ = 180. As the cylindrical particle is not symmetric about

3.3 Computation methods 61

x

y

z

x

y

z

Figure 3.9: The left figure shows a horizontally oriented plate and the rightfigure shows a horizontally oriented cylinder.

the propagation direction, the scattered radiation is polarized, evenfor forward and backward scattering directions. For oriented particlesthe phase matrix does not depend only on the scattering angle buton the incident and scattered directions with respect to the particleorientation. For different directions, the phase matrix for the platesalso deviates strongly from the Mie phase matrix.

This very short analysis of these special phase matrices shows thatparticle shape has a significant impact on the intensity and the polar-ization signal of microwave radiation in the atmosphere if the cirruscloud particles are oriented.

3.3.4 Further methods

Other methods to calculate single scattering properties are the Dis-crete Dipole Approximation (DDA), the geometrical optics approxi-mation, different extended T-matrix methods and several other meth-ods, which are not discussed here.

Cloud particles are not so large that their optical properties in themicrowave region could be calculated in the geometrical optics ap-proximation. Therefore we have not considered this method.

The public domain DDA program DDSCAT developed by Draineand Flateau (2003) is available at http://www.astro.princeton.edu/∼draine/DDSCAT.html. This program is widely used and tested for

62 3 Clouds as scattering media

0 50 100 1500

0.5

1

1.5

2

2.5

Θ [ ° ]

F11

[ ]

Mie aspect ratios: 0.3 0.5 1.0 2.0 5.0

0 50 100 150−1

−0.5

0

0.5

1

1.5

Θ [ ° ]

F12

[ ]

0 50 100 150−0.2

−0.1

0

0.1

0.2

Θ [ ° ]

F34

[ ]

0 50 100 150−2

−1

0

1

2

3

Θ [ ° ]

F33

[ ]

Figure 3.10: Grey line: Phase matrix elements calculated using Mie-theory.Black lines: Phase matrix elements for horizontally oriented plates/cylinderswith different aspect ratios calculated using the T-matrix program for ar-bitrary oriented particles. The dashed lines correspond to cylinders withaspect ratios of 0.3 (thick) and 0.5 (thin). The dotted lines correspond toplates with aspect ratios of 2.0 (thin) and 5.0 (thick). The solid line corre-sponds to a cylinder with an aspect ratio of 1.0. The equal volume sphereradius is 75 µm and the frequency is 318 GHz.

different types of particles. DDA can deal with any shapes of particles.The most noticeable shortcoming of DDA is its tremendous demandin computing time and memory. Compared to DDA, the T-matrixmethod is much faster and has much less demand in computer mem-ory. However the Mishchenko code is designed only for particles withrotational symmetry and therefore it can be used only for a limitednumber of particle types. So far we have not used DDA for ARTSsimulations, but it would be interesting to study the effect of particle

3.4 Single scattering properties ARTS 63

shape for particles which are not rotationally symmetric and/or haveextreme aspect ratios. This is planned in future studies.

The T-matrix results for spherical particles were compared to re-sults obtained using the extended T-matrix code for aggregates, whichwas developed by Havemann and Baran (2001). The result is shown inFigure 3.11. The black lines show scattering, absorption and extinc-tion efficiencies for spherical particles. The circles and cubes show theresults for aggregates with aspect ratios of 1.0 and 4.0 respectively.For frequencies below 400 GHz there are only very small differencesbetween the results of the two different aggregate particle types andthe results for spherical particles. Only for frequencies about 500 GHzthe difference becomes more significant, especially for the aggregatewith an aspect ratio of 4.0. A drawback of the extended T-matrixprogram is that it provides only the cross sections, not the full extinc-tion and phase matrices, which are required for polarized radiativetransfer simulations. For this reason and because of the fact that thecross sections do not deviate much from the cross sections for spheri-cal particles, it was decided not to use the extended T-matrix code asan input for ARTS.

3.4 Single scattering propertiesin the ARTS model

As shown in Section 3.1, clouds consist of a variety of particle sizes andshapes. Furthermore, the cloud particles can be oriented, in most casesthey are horizontally aligned. It is not possible to model the cloudsexactly as they are in nature, therefore we need some approximations.In ARTS different kinds of scattering media are implemented. Onekind consists of randomly oriented particles, which allows very effi-cient computation of single scattering properties. Although this kindof scattering medium is only a special case, it provides a rather goodnumerical description of the scattering properties of clouds and is byfar the most often used theoretical model in particle scattering theory.As the polarization signal of clouds depends significantly on the cloud

64 3 Clouds as scattering media

100 150 200 250 300 350 400 450 500 5500

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

frequency [GHz]

effic

ienc

y [ ]

scattering efficiencyabsorption efficiencyextinction efficiencyaggregate aspect ratio = 1.0aggregate aspect ratio = 4.0

Figure 3.11: Comparison between different T-matrix codes. The black linescorrespond to extinction (solid line), scattering (dashed line) and absorption(dotted line) efficiency obtained for spherical particles using the T-matrixcode for randomly oriented particles. The crosses correspond to results ob-tained for aggregates with an aspect ratio of 1.0 and the circles are theresults for aggregates with an aspect ratio of 4.0.

particle orientation, there is also a kind of scattering media consistingof horizontally aligned particles.

The single scattering properties are pre-calculated using the T-matrix method and stored in data-files. The structure of these data-files is shown in Table 3.1. The first field of the structure includes avalue (enum), which characterizes the kind of scattering medium, i.e.,whether the particles are arbitrarily oriented, randomly oriented orhorizontally aligned. This information is needed, since different scat-tering media types are treated differently inside the program. Thenext field is a string which includes a short description of the scatter-ing medium contained in the data-file. This should include informationabout how the single scattering properties were generated. After that,numerical grids for variables, on which the single scattering propertiesdepend on, are stored. These are the frequency grid, the temperature

3.4 Single scattering properties ARTS 65

grid, the zenith angle grid and the azimuth angle grid. Note that thesame numerical grids are used for incident and scattered directions.The last three fields contain the data. First the phase matrix data isstored as a seven-dimensional array, as the phase matrix depends onthe frequency (ν), the temperature (T ), the scattered direction (θ, φ),the incident direction (θ′, φ′) and it has in general sixteen matrix ele-ments. The extinction matrix data and the absorption coefficient vec-tor data are five-dimensional arrays. They have two dimensions lessthan the array holding the phase matrix data, since they are definedonly for the forward direction.

Table 3.1: Structure of single scattering data files

Symbol Type Dimensions Description

enum specification of particle typeString short description of particle type

ν Vector (ν) frequency gridT Vector (T ) temperature gridθ Vector (θ) zenith angle gridφ Vector (φ) azimuth angle grid

〈Z〉 7D Array (ν, T, θ, φ, θ′, φ′, i) phase matrix〈K〉 5D Array (ν, T, θ, φ, i) extinction matrix〈a〉 5D Array (ν, T, θ, φ, i) absorption vector

The following three kinds of scattering media are implemented inARTS so far. The number in brackets is the corresponding enum-value.General case (p10): If the scattering medium does not have any

symmetries, all sixteen elements of the phase matrix have to bestored. The individual phase matrices are calculated usingMishchenko’s T-matrix code for single particles in fixed orienta-tion, which is described in Section 3.3.3. The extinction matrix hasin general seven independent elements and the absorption vectorhas four different elements.

Randomly oriented particles (p20): For this case the scatteringmedium is macroscopically isotropic and mirror-symmetric and wecalculate the single scattering properties in the particle frame (Fig-

66 3 Clouds as scattering media

ure 3.3) using the T- matrix code for randomly oriented particles,which is described in Section 3.3.3. This reduces the size of thedatafiles enormously, as the single scattering properties depend onlyon the scattering angle, not on four angles needed to describe in-cident and scattered directions. Furthermore the number of inde-pendent elements of the phase matrix, the extinction matrix andthe absorption coefficient vector is reduced. Only six elements ofthe phase matrix are independent and the extinction matrix is di-agonal, therefore only one element needs to be stored in the datafiles. Only the first element of the absorption vector is non-zero.Moreover, extinction and absorption are independent of the prop-agation direction. The only drawback is that the single scatteringdata has to be transformed from the particle frame representationto the laboratory frame representation. These transformations aredescribed in detail in Appendix B.2.

Horizontally aligned plates and columns (p30): For particlesthat are azimuthally randomly oriented, one angular dimension ofthe phase matrix data array is redundant, as the phase matrix isindependent of the incident azimuth angle. Furthermore, regardingthe symmetry of this case, it can be shown that for the scattereddirections only half of the angular grids are required. As for thegeneral case, the fixed orientation T-matrix code for single scat-tering particles is used. The averaging over azimuthal orientationsis done using the exact T -matrix averaging method of Mishchenkoet al. (2000) for the extinction matrix, and by numerical integrationfor the phase matrix. The data is stored in the laboratory frameomitting the redundant data. Therefore this data does not need tobe transformed. In order to use it for the radiative transfer equa-tion, we only need appropriate reading and interpolation routinesfor this data format.It is very convenient to use the PYTHON module PyARTS, which

has been developed especially for ARTS and which is freely available athttp://www.sat.uni-bremen.de/cgi-bin/cvsweb.cgi/PyARTS/. Thismodule can be used to generate single scattering properties for hori-zontally aligned as well as for randomly oriented particles in the ARTSdata-file-format. PyARTS has been developed by C. Davis, who has

3.5 Particle size distributions 67

implemented the Monte Carlo scattering algorithm in ARTS (see Sec-tion 5.3).

3.5 Representation of the particle sizedistribution

The particle size has an important impact on the scattering and ab-sorption properties of cloud particles as shown in Figure 3.6. Cloudscontain a whole range of different particle sizes, which can be describedby a size distribution giving the number of particles per unit volumeper unit radius interval as a function of radius. It is most convenientto parameterize the size distribution by analytical functions, becausein this case optical properties can be calculated much faster than forarbitrary size distributions. The T-matrix code for randomly orientedparticles includes several types of analytical size distributions, e.g.,the gamma distribution or the log-normal distribution. This sectionpresents the size distribution parameterizations, which were used forthe ARTS simulations included in this thesis.

3.5.1 Mono-disperse particle distribution

The most simple assumption is, that all particles in the cloud havethe same size. In order to study scattering effects like polarizationor the influence of particle shape, it makes sense to use this mostsimple assumption, because one can exclude effects coming from theparticle size distribution itself. This simple assumption was made inthe simulations presented in Chapter 7.

Along with the single scattering properties we need the particlenumber density field, which specifies the number of particles per cubicmeter at each grid point, for ARTS scattering simulations. For a givenIMC and mono-disperse particles the particle number density np issimply

np(IMC, r) =IMCm

=IMCV ρ

=34π

IMCρr3

, (3.10)

68 3 Clouds as scattering media

where m is the mass of a particle, r is its equal volume sphere radius,ρ is its density, and V is its volume.

3.5.2 Gamma size distribution

A commonly used distribution for radiative transfer modeling in cirrusclouds is the gamma distribution

n(r) = arα exp(−br). (3.11)

The dimensionless parameter α describes the width of the distribution.The other two parameters can be linked to the effective radius Reff

and the ice mass content IMC as follows:

b =α + 3Reff

, (3.12)

a =IMC

4/3πρb−(α+4)Γ[α + 4], (3.13)

where ρ is the density of the scattering medium and Γ is the gammafunction. For cirrus clouds ρ corresponds to the bulk density of ice,which is 917 kg/m3.

Generally, the effective radius Reff is defined as the average radiusweighted by the particle cross-section

Reff =1〈A〉

∫ rmax

rmin

A(r)rn(r)dr, (3.14)

where A is the area of the geometric projection of a particle. Theminimal and maximal particle sizes in the distribution are given byrmin and rmax respectively. In the case of spherical particles A = πr2.The average area of the geometric projection per particle 〈A〉 is givenby

〈A〉 =

∫ rmax

rminA(r)n(r)dr∫ rmax

rminn(r)dr

. (3.15)

The question is how well a gamma distribution can represent thetrue particle size distribution in radiative transfer calculations. This

3.5 Particle size distributions 69

question is investigated by Evans et al. (1998). The authors cometo the conclusion that a gamma distribution represents the distribu-tion of realistic clouds quite well, provided that the parameters Reff ,IMC and α are chosen correctly. They show that setting α = 1 andcalculating only Reff gave an agreement within 15% in 90% of theconsidered measurements obtained during the First ISCCP RegionalExperiment (FIRE). Therefore, for all calculations including gammasize distributions for ice particles, α = 1 was assumed. The results ofthese calculations are presented in Chapters 5 and 6.

The particle number density for size distributions is obtained byintegration of the distribution function over all sizes:

np(IMC, Reff) =∫ ∞

0

n(r)dr (3.16)

=∫ ∞

0

arα exp(−br)dr = aΓ(α + 1)

bα+1. (3.17)

After setting α = 1, inserting Equation (3.13) and some simple algebrawe obtain

np(IMC, Reff) =2π

IMCρR3

eff

. (3.18)

Comparing Equations (3.10) and (3.18), we see that the particle num-ber density for mono-disperse particles with a particle size of R issmaller than the particle number density for gamma distributed par-ticles with Reff = R. The reason is that in the gamma distributionmost particles are smaller than Reff .

3.5.3 Ice particle size parameterization byMcFarquhar and Heymsfield

A more realistic parameterization of tropical cirrus ice crystal size dis-tributions was derived by McFarquhar and Heymsfield (1997), who de-rived the size distribution as a function of temperature and IMC. Theparameterization was made based on observations during the CentralEquatorial Pacific Experiment (CEPEX). Smaller ice crystals with anequal volume sphere radius of less than 50 µm are parametrized as a

70 3 Clouds as scattering media

sum of first-order gamma functions:

n(r) =12 · IMC<50α

5<50r

πρΓ(5)exp(−2α<50r), (3.19)

where α<50 is a parameter of the distribution, and IMC<50 is the massof all crystals smaller than 50 µm in the observed size distribution.Large ice crystals are represented better by a log-normal function

n(r) =3 · IMC>50

π3/2ρ√

2 exp(3µ>50 + (9/2)σ2>50)rσ>50r3

0

· exp

−12

(log 2r

r0− µ>50

σ>50

)2 , (3.20)

where IMC>50 is the mass of all ice crystals greater than 50 µm in theobserved size distribution, r0 = 1 µm is a parameter used to ensurethat the equation does not depend on the choice of unit for r, σ>50 isthe geometric standard deviation of the distribution, and µ>50 is thelocation of the mode of the log-normal distribution. The fitted param-eters of the distribution can be looked up in the article by McFarquharand Heymsfield (1997). The particle number density field is obtainedby numerical integration over a discrete set of size bins. This pa-rameterization of particle size has been implemented in the PyARTSpackage, which was introduced in Section 3.4. Using PyARTS onecan calculate the size distributions, the corresponding single scatter-ing properties and the particle number density fields for given IMCand temperature. Calculations using this parameterization of cloudparticle sizes are presented in Chapter 8.

4 The DOIT scattering model

The Discrete Ordinate ITerative (DOIT) method is one of the scat-tering algorithms in ARTS. Besides the DOIT method a backwardMonte Carlo scheme has been implemented (see Section 5.3). TheDOIT method is unique because a discrete ordinate iterative methodis used to solve the scattering problem in a 3D spherical atmosphere.A literature review about scattering models for the microwave region,which is presented in Appendix A, shows that former implementa-tions of discrete ordinate schemes are only applicable for (1D-)plane-parallel or 3D-cartesian atmospheres. All of these algorithms can notbe used for the simulation of limb radiances. A description of theDOIT method, similar to what is presented in this chapter, has beenpublished in Emde et al. (2004a).

4.1 The discrete ordinate iterative method4.1.1 Radiation field

The Stokes vector depends on the position in the cloud box and on thepropagation direction specified by the zenith angle (θ) and the azimuthangle (φ). All these dimensions are discretized inside the model; fivenumerical grids are required to represent the radiation field I:

~p = p1, p2, ..., pNp,

~α = α1, α2, ..., αNα,~β = β1, β2, ..., βNβ

, (4.1)~θ = θ1, θ2, ..., θNθ

,~φ = φ1, φ2, ..., φNφ

.

71

72 4 The DOIT scattering model

Here ~p is the pressure grid, ~α is the latitude grid and ~β is the longitudegrid. The radiation field is a set of Stokes vectors (Np × Nα × Nβ ×Nθ ×Nφ elements) for all combinations of positions and directions:

I = I1(p1, α1, β1, θ1, φ1), I2(p2, α1, β1, θ1, φ1), ...,

INp×Nα×Nβ×Nθ×Nφ(pNp

, αNα, βNβ

, θNθ, φNφ

). (4.2)

In the following we will use the notation

i = 1 . . . Np

j = 1 . . . Nα

I = Iijklm k = 1 . . . Nβ . (4.3)

l = 1 . . . Nθ

m = 1 . . . Nφ

4.1.2 Vector radiative transfer equation solution

Figure 4.1 shows a schematic of the iterative method, which is appliedto solve the vector radiative transfer equation (1.42).

The first guess field

I(0) =

I(0)ijklm

, (4.4)

is partly determined by the boundary condition given by the radiationcoming from the clear sky part of the atmosphere traveling into thecloud box. Inside the cloud box an arbitrary field can be chosen as afirst guess. In order to minimize the number of iterations it should beas close as possible to the solution field.

The next step is to solve the scattering integrals⟨S

(0)ijklm

⟩=∫

4π

dn′ 〈Zijklm〉 I(0)ijklm, (4.5)

using the first guess field, which is now stored in a variable reservedfor the old radiation field. For the integration we use equidistant an-gular grids in order to save computation time (cf. Section 4.3). Theradiation field, which is generally defined on finer angular grids (~φ, ~θ),

4.1 The discrete ordinate iterative method 73

field

yes

no

old radiation field

new radiation field solution field

convergence

test

radiative transfer step

(averaged quantities)

scattering intergral

first guess field

Figure 4.1: Schematic of the iterative method to solve the VRTE in thecloud box.

is interpolated on the equidistant angular grids. The integration isperformed over all incident directions n′ for each propagation direc-tion n. The evaluation of the scattering integral is done for all gridpoints inside the cloud box. The obtained integrals are interpolatedon ~φ and ~θ. The result is the first guess scattering integral field S0:

S(0) =⟨

S(0)ijklm

⟩. (4.6)

Figure 4.2 shows a propagation path step from a grid point P =(pi, αj , βk) into direction n = (θl, φm). The radiation arriving at P

from the direction n′ is obtained by solving the linear differentialequation:

dI(1)

ds= −〈K〉I(1) + 〈a〉 B +

⟨S(0)

⟩, (4.7)

where 〈K〉, 〈a〉, B and⟨S(0)

⟩are averaged quantities. This equa-

tion can be solved analytically for constant coefficients. Multi-linearinterpolation gives the quantities K ′,a′,S′ and T ′ at the intersectionpoint P ′. To calculate the radiative transfer from P ′ towards P all

74 4 The DOIT scattering model

( θ , φm )l

( θ , φm )’l’

( j , βk )pi , α

, j−1α , β k)i−1p(

( j , βk )pi , α’ ’ ’

p( i+1, α ,β k)j p( i+1, α j+1, β

p( i α j+1, β k,

k

)

)

Figure 4.2: Path from a grid point ((pi, αj , βk) - (×)) to the intersectionpoint ((p′i, α′

j , β′k) - ()) with the next grid cell boundary. Viewing direction

is specified by (θl, φm) at (×) or (θ′l, φ′m) at ().

quantities are approximated by taking the averages between the val-ues at P ′ and P . The average value of the temperature is used to getthe averaged Planck function B.

The solution of Equation (4.7) is found analytically using a matrixexponential approach (see Appendix B.1):

I(1) = e−〈K〉sI(0)(I− e−〈K〉s

)〈K〉

−1(〈a〉 B +

⟨S(0)

⟩), (4.8)

where I denotes the identity matrix and I(0) the initial Stokes vector.The radiative transfer step from P ′ to P is calculated, therefore I(0)

is the incoming radiation at P ′ into direction (θ′l, φ′m), which is the

first guess field interpolated on P ′. This radiative transfer step calcu-lation is done for all points inside the cloud box in all directions. Theresulting set of Stokes vectors (I(1) for all points in all directions) isthe first order iteration field I(1):

I(1) =

I(1)ijklm

. (4.9)

The first order iteration field is stored in a variable reserved for thenew radiation field.

4.1 The discrete ordinate iterative method 75

In the convergence test the new radiation field is compared to the oldradiation field. For the difference field, the absolute values of all Stokesvector elements for all cloud box positions are calculated. If one of thedifferences is larger than a requested accuracy limit, the convergencetest is not fulfilled. The user can define different convergence limitsfor the different Stokes components.

If the convergence test is not fulfilled, the first order iteration fieldis copied to the variable holding the old radiation field, and is thenused to evaluate again the scattering integral at all cloud box points:⟨

S(1)ijklm

⟩=∫

4π

dn′ 〈Z〉 I(1)ijklm. (4.10)

The second order iteration field

I(2) =

I(2)ijklm

, (4.11)

is obtained by solving

dI(2)

ds= −〈K〉I(2) + 〈a〉 B +

⟨S(1)

⟩, (4.12)

for all cloud box points in all directions. This equation contains alreadythe averaged values and is valid for specified positions and directions.

As long as the convergence test is not fulfilled the scattering integralfields and higher order iteration fields are calculated alternately.

We can formulate a differential equation for the n-th order iterationfield. The scattering integrals are given by⟨

S(n−1)ijklm

⟩=∫

4π

dn′ 〈Z〉 I(n−1)ijklm , (4.13)

and the differential equation for a specified grid point into a specifieddirection is

dI(n)

ds= −〈K〉I(n) + 〈a〉 B +

⟨S(n−1)

⟩. (4.14)

Thus the n-th order iteration field

I(n) =

I(n)ijklm

, (4.15)

76 4 The DOIT scattering model

is given by

I(n) = e−〈K〉s · I(n−1)(I− e−〈K〉s)〈K〉−1

(〈a〉 B +⟨S(n−1)

⟩),

(4.16)

for all cloud box points and all directions defined in the numericalgrids.

If the convergence test∣∣∣I(N)ijklm (pi, αj , βk, θl, φm)− I

(N−1)ijklm (pi, αj , βk, θl, φm)

∣∣∣ < ε,

(4.17)

is fulfilled, a solution to the vector radiative transfer equation (1.42)has been found:

I(N) =

I(N)ijklm

. (4.18)

4.1.3 Scalar radiative transfer equation solution

In analogy to the scattering integral vector field the scalar scatteringintegral field is obtained:⟨

S(0)ijklm

⟩=∫

4π

dn′ 〈Z11〉 I(0)ijklm. (4.19)

The scalar radiative transfer equation (1.51) with a fixed scatteringintegral is

dI(1)

ds= −〈K11〉 I(1) + 〈a1〉B +

⟨S(0)

⟩. (4.20)

Assuming constant coefficients this equation is solved analytically af-ter averaging extinction coefficients, absorption coefficients, scatteringvectors and the temperature. The averaging procedure is done analo-gously to the procedure described for solving the VRTE. The solutionof the averaged differential equation is

I(1) = I(0)e−〈K11〉s +〈a1〉 B +

⟨S(0)

⟩〈K11〉

(1− e−〈K11〉s

), (4.21)

4.1 The discrete ordinate iterative method 77

where I(0) is obtained by interpolating the initial field, and 〈K11〉, 〈a1〉,B and

⟨S(0)

⟩are the averaged values for the extinction coefficient, the

absorption coefficient, the Planck function and the scattering integralrespectively. Applying this equation leads to the first iteration scalarintensity field, consisting of the intensities I(1) at all points in thecloud box for all directions.

As the solution to the vector radiative transfer equation, the solu-tion to the scalar radiative transfer equation is found numerically bythe same iterative method. The convergence test for the scalar equa-tion compares the values of the calculated intensities of two successiveiteration fields.

4.1.4 Single scattering approximation

The DOIT method uses the single scattering approximation, whichmeans that for one propagation path step the optical depth is assumedto be much less than one so that multiple-scattering can be neglectedalong this propagation path step. It is possible to choose a rathercoarse grid inside the cloud box. The user can define a limit for themaximum propagation path step length. If a propagation path stepfrom one grid cell to the intersection point with the next grid cellboundary is greater than this value, the path step is divided in severalsteps such that all steps are less than the maximum value. The userhas to make sure that the optical depth due to cloud particles forone propagation path sub-step is is sufficiently small to assume singlescattering. The maximum optical depth due to ice particles is

τmax = 〈Kp〉 ·∆s, (4.22)

where ∆s is the length of a propagation path step. In all simulationspresented in this thesis τmax 0.01 is assumed. This threshold valueis also used in Czekala (1999a). The radiative transfer calculation isdone along the propagation path through one grid cell. All coefficientsof the VRTE are interpolated linearly on the propagation path points.

78 4 The DOIT scattering model

4.2 Sequential updateIn the previous Section 4.1 the iterative solution method for the VRTEhas been described. For each grid point inside the cloud box the inter-section point with the next grid cell boundary is determined in eachviewing direction. After that, all the quantities involved in the VRTEare interpolated onto this intersection point. As described in the sec-tions above, the intensity field of the previous iteration is taken toobtain the Stokes vector at the intersection point. Suppose that thereare N pressure levels inside the cloud box. If the radiation field isupdated taking into account for each grid point only the adjacentgrid cells, at least N -1 iterations are required until the scattering ef-fect from the lower-most pressure level has propagated throughoutthe cloud box up to the uppermost pressure level. From these con-siderations, it follows, that the number of iterations depends on thenumber of grid points inside the cloud box. This means that the orig-inal method is very ineffective where a fine resolution inside the cloudbox is required to resolve the cloud inhomogeneities.

A solution to this problem is the “sequential update of the radiationfield”, which is shown schematically in Figure 4.3. For simplicity it willbe explained in detail for a 1D cloud box. We divide the update ofthe radiation field, i.e., the radiative transfer step calculations for allpositions and directions inside the cloud box, into three parts: Updatefor “up-looking” zenith angles (0 ≤ θup ≤ 90), for “down-looking”angles (θlimit ≤ θdown ≤ 180) and for “limb-looking” angles (90 <

θlimb < θlimit). The “limb-looking” case is needed, because for anglesbetween 90 and θlimit the intersection point is at the same pressurelevel as the observation point. The limiting angle θlimit is calculatedgeometrically. Note that the propagation direction of the radiation isopposite to the viewing direction or the direction of the line of sight,which is indicated by the arrows. In the 1D case the radiation field isa set of Stokes vectors each of which depend upon the position anddirection:

I = I (pi, θl) . (4.23)

The boundary condition for the calculation is the incoming radiation

4.2 Sequential update 79

limb

limb

limb

limb

down

down

up

up

up

θ

θ

θ

θ

θ

θ

θ

θθ

pN

p0

θdown

Figure 4.3: Schematic of the sequential update (1D) showing the three dif-ferent parts: “up-looking” corresponds to zenith angles θup, “limb-looking”corresponds to θlimb “down-looking” corresponds to θdown.

field on the cloud box boundary Ibd:

Ibd = I (pi, θl) where pi = pN ∀ θl ∈ [0, θlimit]

pi = p0 ∀ θl ∈ (θlimit, 180], (4.24)

where p0 and pN are the pressure coordinates of the lower and uppercloud box boundaries respectively. For down-looking directions, theintensity field at the lower-most cloud box boundary and for up- andlimb-looking directions the intensity field at the uppermost cloud boxboundary are the required boundary conditions respectively.

4.2.1 Up-looking directions

The first step of the sequential update is to calculate the intensityfield for the pressure coordinate pN−1, the pressure level below theuppermost boundary, for all up-looking directions. Radiative transfersteps are calculated for paths starting at the uppermost boundary andpropagating to the (N − 1) pressure level. The required input for thisradiative transfer step are the averaged coefficients of the uppermostcloud box layer and the Stokes vectors at the uppermost boundaryfor all up-looking directions. These are obtained by interpolating the

80 4 The DOIT scattering model

boundary condition Ibd on the appropriate zenith angles. Note thatthe zenith angle of the propagation path for the observing directionθl does not equal θ′l at the intersection point due to the sphericalgeometry. If θl is close to 90 this difference is most significant.

To calculate the intensity field for the pressure coordinate pN−2, werepeat the calculation above. We have to calculate a radiative transferstep from the (N − 1) to the (N − 2) pressure level. As input we needthe interpolated intensity field at the (N − 1) pressure level, whichhas been calculated in the last step.

For each pressure level (m − 1) we take the interpolated field ofthe layer above (I(pm)(1)). Using this method, the scattering influ-ence from particles in the upper-most cloud box layer can propagateduring one iteration down to the lower-most layer. This means thatthe number of iterations does not scale with the number of pressurelevels, which would be the case without sequential update.

The radiation field at a specific point in the cloud box is obtainedby solving Equation (4.8). For up-looking directions at position pm−1

we may write:

I (pm−1, θup)(1) = e−〈K(θup)〉sI (pm, θup)(1)

+(I− e−〈K(θup)〉s

)〈K(θup)〉

−1(〈a(θup)〉 B +

⟨S (θup)(0)

⟩).

(4.25)

For simplification we write

I(pm−1, θup)(1) = A(θup)I (pm, θup)(1) + B(θup). (4.26)

Solving this equation sequentially, starting at the top of the cloudand finishing at the bottom, we get the updated radiation field for allup-looking angles.

I(pi, θup)(1) =

I(1) (pi, θl)

∀ θl ∈ [0, 90]. (4.27)

4.2.2 Down-looking directions

The same procedure is done for down-looking directions. The onlydifference is that the starting point is the lower-most pressure level

4.2 Sequential update 81

p1 and the incoming clear sky field at the lower cloud box boundary,which is interpolated on the required zenith angles, is taken as bound-ary condition. The following equation is solved sequentially, startingat the bottom of the cloud box and finishing at the top:

I(pm, θdown)(1) = A(θdown)I (pm−1, θdown)(1) + B(θdown). (4.28)

This yields the updated radiation field for all down-looking angles.

I(pi, θdown)(1) =

I(1) (pi, θl)

∀ θl ∈ [θlimit, 180]. (4.29)

4.2.3 Limb directions

A special case for limb directions, which correspond to angles slightlyabove 90 had to be implemented. If the tangent point is part ofthe propagation path step, the intersection point is exactly at thesame pressure level as the starting point. In this case the linearlyinterpolated clear sky field is taken as input for the radiative transfercalculation, because we do not have an already updated field for thispressure level:

I(pm, θlimb)(1) = A(θlimb)I (pm, θlimb)(0) + B(θlimb) (4.30)

By solving this equation the missing part of the updated radiationfield is obtained

I(pi, θlimb)(1) = I (pi, θl) ∀ θl ∈]90, θlimit[ (4.31)

For all iterations the sequential update is applied. Using this methodthe number of iterations depends only on the optical thickness of thecloud or on the number of multiple-scattering events, not on the num-ber of pressure levels.

How the sequential update is performed in the 3D model is describedin Eriksson et al. (2004).

82 4 The DOIT scattering model

4.3 Grid optimization and interpolationmethods

The accuracy of the DOIT method depends very much on the dis-cretization of the zenith angle. The reason is that the intensity fieldstrongly increases at about θ = 90. For angles below 90 (“up-looking” directions) the intensity is very small compared to anglesabove 90 (“down-looking” directions), because the thermal emis-sion from the lower atmosphere and from the ground is much largerthan thermal emission from trace gases in the upper atmosphere. Fig-ure 4.4 shows an example intensity field as a function of zenith anglefor different pressure levels inside a cloud box, which is placed from7.3 to 12.7 km altitude, corresponding to pressure limits of 411 hPaand 188 hPa respectively. The cloud box includes 27 pressure levels.The frequency of the sample calculation was 318 GHz. A midlatitude-summer scenario including water vapor, ozone, nitrogen and oxygenwas used. The atmospheric data was taken from the FASCOD (An-derson et al., 1986) and the spectroscopic data was obtained fromthe HITRAN database (Rothman et al., 1998). For simplicity this 1Dset-up was chosen for all sample calculations in this section. As theintensity (or the Stokes vector) at the intersection point of a prop-agation path is obtained by interpolation, large interpolation errorscan occur for zenith angles of about 90 if the zenith angle grid dis-cretization is too coarse. Taking a very fine equidistant zenith anglegrid leads to very long computation times. Therefore a zenith anglegrid optimization method is required.

For the computation of the scattering integral it is possible to take amuch coarser zenith angle resolution without losing accuracy. It doesnot make sense to use the zenith angle grid, which is optimized torepresent the radiation field with a certain accuracy. The integrandis the product of the phase matrix and the radiation field. The peaksof the phase matrices can be at any zenith angle, depending on theincoming and the scattered directions. The multiplication smoothsout both the radiation field increase at 90 and the peaks of the phasematrices. Test calculations have shown that an increment of 10 issufficient. Taking the equidistant grid saves the computation time of

4.3 Grid optimization and interpolation methods 83

the scattering integral to a very large extent, because much less gridpoints are required.

0 20 40 60 80 100 120 140 160 1800

50

100

150

200

250

300

Zenith angle [ ° ]

BT

[ K ]

p = 411 hPap = 346 hPap = 290 hPap = 242 hPa

Figure 4.4: Intensity field for different pressure levels.

4.3.1 Zenith angle grid optimization

As a reference field for the grid optimization the DOIT method is ap-plied for an empty cloud box using a very fine zenith angle grid. Thegrid optimization routine finds a reduced zenith angle grid which canrepresent the intensity field with the desired accuracy. It first takes theradiation at 0 and 180 and interpolates between these two pointson all grid points contained in the fine zenith angle grid for all pres-sure levels. Then the differences between the reference radiation fieldand the interpolated field are calculated. The zenith angle grid point,where the difference is maximal is added to 0 and 180. After thatthe radiation field is interpolated between these three points formingpart of the reduced grid and again the grid point with the maximumdifference is added. Using this method more and more grid points are

84 4 The DOIT scattering model

0 20 40 60 80 100 120 140 160 1800

100

200

300

Zenith angle [ ° ]

BT

[ K ]

clearsky

90 95 100 105−0.5

0

0.5

1

Zenith angle [ ° ]

Inte

rpol

atio

n er

ror [

% ]

acc: 0.1%acc: 0.2%acc: 0.5%

100 120 140 160 180−0.15

−0.1

−0.05

0

0.05

0.1

Zenith angle [ ° ]

Inte

rpol

atio

n er

ror [

% ]

acc: 0.1%acc: 0.2%acc: 0.5%

Figure 4.5: Interpolation errors for different grid accuracies. Top panel:Clear sky radiation simulated for a sensor at an altitude of 13 km forall viewing directions. Bottom left: Grid optimization accuracy for limbdirections. Bottom right: Grid optimization accuracy for down-lookingdirections.

added to the reduced grid until the maximum difference is below arequested accuracy limit.

The top panel of Figure 4.5 shows the clear sky radiation in allviewing directions for a sensor located at 13 km altitude. This resultwas obtained with a switched-off cloud box. The difference betweenthe clear sky part of the ARTS model and the scattering part is thatin the clear sky part the radiative transfer calculations are done alongthe line of sight of the instrument whereas inside the cloud box the RTcalculations are done as described in the previous section to obtainthe full radiation field inside the cloud box. In the clear sky part theradiation field is not interpolated, therefore we can take the clear skysolution as the exact solution.

4.3 Grid optimization and interpolation methods 85

The interpolation error is the relative difference between the exactclear sky calculation (cloud box switched off) and the clear sky calcu-lation with empty cloud box. The bottom panels of Figure 4.5 showthe interpolation errors for zenith angle grids optimized with threedifferent accuracy limits (0.1%, 0.2% and 0.5%.). The left plot showsthe critical region close to 90. For a grid optimization accuracy of0.5% the interpolation error becomes very large, the maximum erroris about 8%. For grid accuracies of 0.2% and 0.1% the maximum in-terpolation errors are about 0.4% and 0.2% respectively. However formost angles it is below 0.2%, for all three cases. For down-lookingdirections from 100 to 180 the interpolation error is at most 0.14%for grid accuracies of 0.2% and 0.5% and for a grid accuracy of 0.1%it is below 0.02%.

4.3.2 Interpolation methods

Two different interpolation methods can be chosen in ARTS for theinterpolation of the radiation field in the zenith angle dimension: lin-ear interpolation or three-point polynomial interpolation. The poly-nomial interpolation method produces more accurate results providedthat the zenith angle grid is optimized appropriately. The linear in-terpolation method on the other hand is safer. If the zenith anglegrid is not optimized for polynomial interpolation one should use thesimpler linear interpolation method. Apart from the interpolation ofthe radiation field in the zenith angle dimension linear interpolationis used everywhere in the model. Figure 4.6 shows the interpolationerrors for the different interpolation methods. Both calculations areperformed on optimized zenith angle grids, for polynomial interpo-lation 65 grid points were required to achieve an accuracy of 0.1%and for linear interpolation 101 points were necessary to achieve thesame accuracy. In the region of about 90 the interpolation errors arebelow 1.2% for linear interpolation and below 0.2% for polynomial in-terpolation. For the other down-looking directions the differences arebelow 0.08% for linear and below 0.02% for polynomial interpolation.It is obvious that polynomial interpolation gives more accurate results.

86 4 The DOIT scattering model

Another advantage is that the calculation is faster because less gridpoints are required, although the polynomial interpolation method it-self is slower than the linear interpolation method. Nevertheless, wehave implemented the polynomial interpolation method so far only inthe 1D model. In the 3D model, the grid optimization needs to bedone over the whole cloud box, where it is not obvious that one cansave grid points. Applying the polynomial interpolation method usingnon-optimized grids can yield much larger interpolation errors thanthe linear interpolation method.

90 95 100 105−1.2

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

Zenith angle [ ° ]

Inte

rpol

atio

n Er

ror [

% ]

linearpolynomial

100 120 140 160 180−0.04

−0.02

0

0.02

0.04

0.06

0.08

Zenith angle [ ° ]

Inte

rpol

atio

n Er

ror [

% ]

linearpolynomial

Figure 4.6: Interpolation errors for polynomial and linear interpolation.

4.3.3 Error estimates

The interpolation error for scattering calculations can be estimatedby comparison of a scattering calculation performed on a very finezenith angle grid (resolution 0.001 from 80 to 100) with a scat-tering calculation performed on an optimized zenith angle grid with0.1% accuracy. The cloud box used in previous test calculations isfilled with spheroidal particles with an aspect ratio of 0.5 from 10 to12 km altitude. The ice mass content is assumed to be 4.3 · 10−3 g/m3

at all pressure levels. An equal volume sphere radius of 75 µm is as-sumed. The particles are either completely randomly oriented (p20)or azimuthally randomly oriented (p30) (cf. Appendix 3.4). The top

4.3 Grid optimization and interpolation methods 87

panels of Figure 4.7 show the interpolation errors of the intensity. Forboth particle orientations the interpolation error is in the same rangeas the error for the clear sky calculation, below 0.2 K. The bottompanels show the interpolation errors for Q. For the randomly orientedparticles the error is below 0.5%. For the horizontally aligned parti-cles with random azimuthal orientation it goes up to 2.5% for a zenithangle of about 91.5. It is obvious that the interpolation error for Q

must be larger than that for I because the grid optimization is ac-complished using only the clear-sky field, where the polarization iszero. Only the limb directions about 90 are problematic, for otherdown-looking directions the interpolation error is below 0.2%.

88 4 The DOIT scattering model

90 95 100 105−0.2

−0.1

0

0.1

0.2

0.3

Zenith angle [ ° ]

Inte

rpol

atio

n E

rror

I [ %

]

p20p30

100 120 140 160 180−0.04

−0.03

−0.02

−0.01

0

0.01

0.02

0.03

Zenith angle [ ° ]

Inte

rpol

atio

n E

rror

I [ %

]

p20p30

90 95 100 105−1.5

−1

−0.5

0

0.5

1

1.5

2

2.5

Zenith angle [ ° ]

Inte

rpol

atio

n E

rror

Q [

% ]

p20p30

100 120 140 160 180−0.2

−0.15

−0.1

−0.05

0

0.05

0.1

0.15

0.2

Zenith angle [ ° ]

Inte

rpol

atio

n E

rror

Q [

% ]

p20p30

Figure 4.7: Interpolation errors for a scattering calculation. Left panels: In-terpolation errors for limb directions. Right panels: Interpolation errors fordown-looking directions. Top: Intensity I, Bottom: Polarization difference Q

5 Comparison of the DOIT methodwith other scattering models

5.1 Comparison with the pseudo-sphericalunpolarized model FM2D

An intercomparison between the two-dimensional forward modelFM2D developed at RAL (Rutherford Appleton Laboratory) and theARTS-DOIT scattering model was performed in order to validate theFM2D model, which includes several approximations for efficiencyreasons. The focus in the FM2D development was, that the modelshould be sufficiently fast for use in non-linear retrieval simulations.

5.1.1 The pseudo-spherical approach

The RAL forward model is briefly described and validated in von Clar-mann et al. (2003). The model allows modeling of the atmosphere intwo dimensions with the field expressed as a two-dimensional, verticaland horizontal, section. Operation in 1D mode, where the atmosphereis assumed to be spherically symmetric, is also possible. The modelcan be used for radiative transfer calculations in the microwave and inthe mid-infrared wavelength regions. Scattering has been included byusing the plane-parallel version of the GOMETRAN model (Rozanovet al., 1997), which was extended to include thermal emission as wellas solar radiation. This allows to efficiently calculate the multiple-scattering source function at all points along the limb line-of-sight,ray-traced through a spherical atmosphere. The implementation ofthe scattering in FM2D is described in Kerridge et al. (2004), wherethe following intercomparison study has also been published.

89

90 5 Scattering model intercomparisons

Integration of the RTE

As the full solution of the VRTE (1.42) is computationally expensive,the following assumptions are included in FM2D for simplification:1. Polarization can be neglected, i.e., the SRTE (1.51) is considered.2. The scattering integral can be evaluated using a plane-parallel, 1D-

model with optical properties taken from the 2D-atmosphere at thehorizontal position of the considered propagation path point. Hav-ing obtained the scattering integral, which is often called sourcefunction, at all points along the line-of-sight, the scalar radia-tive transfer equation with a fixed scattering integral term Equa-tion (4.20) is evaluated along the LOS.

The scattering integral is modeled using a modified version of theGOMETRAN++ forward model. The integration of the radiativetransfer equation is done numerically by iterating the following equa-tion, which follows from Equation (1.51), along the line of sight:

Il+1 = Ile−〈K11,l〉∆sl +

〈a1, l〉 B + 〈Sl〉〈K11,l〉

(1− e−〈K11,l〉∆sl

), (5.1)

where l is the index of a path segment and ∆sl the path length. Thecalculation of the source function 〈S〉 is done for all along track gridpoints as a first step.

⟨Sl

⟩is then obtained by interpolating 〈S〉 in

altitude, horizontal along-track position and local LOS zenith angle.

Representation of scatterers

The distribution of scatterers in FM2D is defined by the user as the 2Ddistributions of the associated scattering and absorption coefficient (inkm−1), together with parameters describing the phase function. Thescattering properties are calculated outside FM2D using for examplea Mie scattering program.

In principle the phase function must be specified in FM2D as a func-tion of altitude, along-track position, frequency and scattering angle.In order to minimize the number of input parameters, the angulardependence of the phase function is modeled by two parameters: theHenry-Greenstein asymmetry parameter and the Rayleigh fraction.

5.1 FM2D - A 2D pseudo-spherical model 91

Notes on accuracy of the model

For 1D geometries, the plane-parallel approach is well known to be ac-curate for near-nadir geometries. The extension to limb geometry inFM2D should give similar accurate results, since the plane parallel ap-proximation is only made to determine the scattering integral, wherean integration is performed over all spatial directions. Therefore theaccuracy of the scattering integral will be similar for nadir and limbviewing geometry. This argument does not hold for strongly peakedphase functions which imply strong forward scattering. In this case,near-limb incident directions dominate the scattering integral and arenot well modeled in the plane-parallel approach.

For 2D calculations many important aspects of the radiative transferalong the line of sight are captured by the approach since in mostcases the local vertical atmospheric state profiles are dominant in thescattering integral calculation.

5.1.2 Clear sky comparison

Before comparing the scattering models it had to be assured that theclear sky calculations agree within a required accuracy.

In these calculations the same atmospheric profiles, gas absorptioncoefficients and general model settings were used. The MASTER-Cband was chosen for the intercomparison as it includes different levelsof gas absorption, a low gas absorption region at about 318 GHz andstrong gas absorption at about 324 GHz.

The comparison is based on the standard profiles used at RAL.The species in the atmosphere are restricted to H2O and O3 and theearth is considered to be spherical with a radius of 6367.62 km. Allatmospheric fields are defined on a 1 km vertical grid, irrespectiveof forward model internal grids, and extend to a height of 50 km.The assumed satellite altitude is 820 km. We have chosen a coarseoptimized frequency grid including 70 frequencies. This is sufficientfor an intercomparison of the scattering models, since the scatteringproperties do not change much in the MASTER band.

The results of the intercomparison, which are presented in Fig-

92 5 Scattering model intercomparisons

310320330−10

1

TH = 1 km TH = 8 km TH = 10 km TH = 11.5 km

318 320 322 324

140

180

220

clear sky BT [K]

Frequency [ GHz]318 320 322 324

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

0.6Rel. difference [%]

Frequency [ GHz]

ARTS (black) vs. FM2D (grey) clear sky

Figure 5.1: Clear sky radiative transfer comparison at different tangentheights (TH).

ure 5.1, show that the level of agreement between the two modelsis very good. For a propagation path step length of 1 km, the high-est relative difference is about 0.6% for a non-refractive atmosphere,this corresponds to a difference of approximately 1 K in brightnesstemperature.

5.1.3 Comparison for cloudy scenarios

Definition of cloud scenarios

The cloud setup for the calculations is summarized in Table 5.1. Thetable includes the ice mass content IMC, the effective radius Reff ofthe size distributions, and the altitude of the different cloud scenarios.The cloud particle size was parametrized according to the gammadistribution introduced in Section 3.5.2. All particles were assumed

5.1 FM2D - A 2D pseudo-spherical model 93

to be spherical. According to FIRE measurements, these scenarioscorrespond to realistic cloud cases.

Table 5.1: Definition of cloud scenarios for the scattering model intercom-parison (DOIT versus FM2D)

Cloud IMC [g m−3] Reff [µm] Altitude [km]

1 0.0001 21.5 10 – 12.12 0.004 34.0 10 – 12.13 0.02 68.5 10 – 12.1

4 0.04 85.5 10 – 12.15 0.1 128.5 10 – 12.1

Results and discussion

The result of the calculations is presented in Figure 5.2. The rowsin the figure correspond to the five scenarios. The left column showsthe simulated radiances at tangent heights of 1 km, 8 km, 10 km and11.5 km. The black lines corresponds to the ARTS-DOIT results andthe grey lines to the FM2D results. The middle column shows thescattering effect, which is the difference between cloudy radiances andclear sky radiances. A positive value corresponds to a brightness tem-perature enhancement due to cloud and a negative value correspondsto a brightness temperature depression. Again, the black lines are theARTS results and the grey lines the FM2D results. The plots show,that the modeled radiances are very similar. The right column in thefigure shows the differences between ARTS and FM2D, more precisely,the difference between the simulated scattering effects.

For scenario 1, the weakest cloud scenario, the difference betweenthe models is below 0.004 K. The thicker the cloud the larger thecloud signal and also the difference between the models. The maximumdifference for scenario 2 is approximately 0.3 K. For scenario 3 and4 the maximum difference is about 1.6 K and for scenario 5 about3.6 K. In scenarios 3 and 4 the FM2D results show larger brightnesstemperatures whereas in scenario 5 the ARTS brightness temperaturesare higher. Figure 5.3 shows the normalized phase functions for the

94 5 Scattering model intercomparisons

310320330−50

5x 10−3

TH = 1 km TH = 8 km TH = 10 km TH = 11.5 km

100

140

180

220

scen

ario

1

cloudy BT [K]

−0.05

0

0.05

0.1

0.15∆ BT [K]

−0.003

−0.001

0.001

0.003

Difference [K]

100

140

180

220

scen

ario

2

−4

−1

2

5

8

−0.3

−0.2

−0.1

0

0.1

140

180

220

scen

ario

3

−50

−30

−10

10

30

−1.1

−0.8

−0.5

−0.2

0.1

100

140

180

220

scen

ario

4

−90

−50

−10

30

−1.3

−0.8

−0.3

0.2

318 320 322 324100

140

180

220

scen

ario

5

Frequency [ GHz]318 320 322 324

−140

−110

−80

−50

−20

Frequency [ GHz]318 320 322 324

0.30.81.31.82.32.83.3

Frequency [ GHz]

ARTS (black) vs. FM2D (grey)

Figure 5.2: Results for all scenarios. Left: Cloudy radiances. Centre: Differ-ences between cloudy and clear sky brightness temperatures. Right: Differ-ences between ARTS and FM2D.

5.1 FM2D - A 2D pseudo-spherical model 95

0 20 40 60 80 100 120 140 160 1800

0.5

1

1.5

2

2.5

3

3.5

4

scattering angle [°]

norm

aliz

ed p

hase

func

tion

[ ]

scenario 1scenario 2scenario 3scenario 4scenario 5

Figure 5.3: Normalized phase functions for the simulated cloud scenarios

modeled cloud scenarios. For scenario 1 and 2, the phase functionshows a forward (0) and a backward (180) peak. The shape of thephase function is rather flat, that means radiation is scattered intoall directions. For scenarios 3 to 5 most of the radiation is scatteredinto the forward directions. In scenario 5 the forward peak is mostpronounced. The more radiation is scattered into forward direction,the less accurate is the plane parallel approximation for the calculationof the scattering integral. Both the increase in cloud optical depth andthe more pronounced forward peak in the phase function explain theincreasing differences between ARTS and RAL-FM2D with increasingice mass content and increasing particle size. Simulations for similarcloud scenarios are also presented in Chapter 6, where the effect ofdifferent cloud parameters is discussed in detail using the ARTS-DOITmodel.

96 5 Scattering model intercomparisons

5.2 Comparison with the single scatteringmodel KOPRA for IR wavelengths

This section presents the intercomparison of the DOIT multiple scat-tering algorithm with the KOPRA (Karlsruhe Optimized and PreciseRadiative Transfer Algorithm) model, which includes a single scatter-ing approach. This intercomparison has been published in Hoepfnerand Emde (2005). The validity of single scattering radiative transfercalculations for simulations of limb emission measurements of cloudsin the mid-infrared spectral region was investigated. This study as-sesses the applicability of relatively fast single scattering calculations,which are important for data analysis of measurements of polar strato-spheric and of cirrus clouds by current and future satellite bornespectrally high-resolution limb-emission sounders. Such instrumentsare for instance MIPAS (Michelson Interferometer for Passive Atmo-spheric Sounding) on Envisat (Fischer and Oelhaf, 1996), launched inMarch 2002, or TES (Troposhperic Emission Spectrometer) on EOS-Aura (Beer et al., 2001), launched in July 2004.

5.2.1 Zero- and single scattering solutions

KOPRA was especially developed for the analysis of spectrally highresolved remote sensing measurements of the earth’s atmosphere inthe mid-infrared (Stiller, 2000; Stiller et al., 2002). The part of themodel describing radiative transfer in the gaseous atmosphere hasbeen validated extensively (Glatthor et al., 1999; von Clarmann et al.,2002, 2003; Tjemkes et al., 2003). Based on a layer-by-layer approachKOPRA models the radiative transfer as a succession of extinction,emission and scattering in homogeneous layers.

From the analytic solution of Equation (1.51) for a homogeneouslayer of thickness s with fixed scattering integral (4.21) we can derivethe discretized radiative transfer equation. If the instrumental line-of-sight traverses L layers it reads:

5.2 KOPRA - A single scattering model for the IR 97

I(sobs) = I(0)L∏

l=1

τ(l) +L∑

l=1

〈a1,l〉Bl + 〈Sl〉〈K11,l〉

(1− τl)L∏

j=l+1

τj

,

(5.2)

where l the index of the layer and τl = exp(−〈K11,l〉 sl) is the trans-mission of layer l with thickness sl. For the determination of the scat-tering integral Sl =

∫4π

dn′〈Z11,l(n,n′)〉Il(n′), the incoming radi-ances are calculated neglecting the scattering source term:

I(s′) = I(0)L′∏

l′=1

τ(l′) +L′∑

l′=1

〈a1,l′〉Bl′

〈K11,l′〉(1− τl′)

L′∏j=l′+1

τj

. (5.3)

The prime symbol denotes that the variables belong to the first orderscattering rays.

Below, four different options of scattering in KOPRA will be com-pared with the ARTS DOIT scattering algorithm:

KOPRA(0) Zero scattering scheme neglecting the scattering sourceterm 〈Sl〉 in Equation (5.2).

KOPRA(1) Zero scattering scheme neglecting the scattering sourceterm 〈Sl(n)〉 and replacing the absorption coefficient 〈a1,l〉 by theextinction coefficient 〈K11,l〉 in Equation (5.2). This increases thethermal source term and might compensate for omitting the scat-tering source term. For optically thick clouds this approach shouldresult in the blackbody radiation emitted from the top of the cloud,since extinction and emission are equal.

KOPRA(2) Single scattering scheme using Equation (5.2) and Equa-tion (5.3).

KOPRA(3) Single scattering scheme using Equation (5.2) and Equa-tion (5.3) in which 〈a1,l′〉 is replaced by 〈K11,l′〉, which could possi-bly compensate for neglecting the multiple scattered component ofthe radiation field.

98 5 Scattering model intercomparisons

5.2.2 Definition of scenarios

In order to avoid any problems with differences in line-by-line absorp-tion calculations, absorption cross-sections were calculated using theKOPRA model. The spectral interval 946.149–950.837 cm−1 was se-lected and the gases CO2 and H2O were taken into account. The cross-sections were calculated for the atmospheric pressure-temperatureprofile on a 0.5 km grid. These pre-calculated cross-sections were usedto simulate the gaseous radiance contributions with the ARTS model.

The altitude profile of the cloud was defined between 9.5 and12.5 km altitude with linearly increasing (from 0 cm−3) values of theparticle number density from 9.5 to 10 km and decreasing values (to0 cm−3) from 12 to 12.5 km. Between 10 and 12 km the number densitywas constant. A log-normal size distribution of spherical ice-cloud par-ticles was assumed with a median-radius of 4 µm and a width of 0.3.Table 5.2 summarizes the five scenarios of increasing density. Here,the smallest volume density is in the order of that of typical polarstratospheric clouds of type I containing a large fraction of HNO3 andscenario 2 is representative of polar stratospheric clouds of type IIconsisting of ice particles.

For the determination of single scattering properties, Mie calcula-tions were done on basis of refractive indexes of ice by Toon et al.(1994). In the middle of the defined spectral interval the refractiveindex is (1.07+0.17i) leading to an absorption cross-section of 5.1 ×10−7 cm2 and a scattering cross-section of 1.6×10−7 cm2. This resultsin a single scattering albedo ω0 = 0.24. Hence this is a case of relativelystrong absorption, since the chosen wave-number region is situated atthe edge of an ice absorption peak in the mid-IR with a maximumaround 830 cm−1. To cover also a case of strong scattering and weakabsorption we used an index of refraction of (1.25+0.018i) for thesame wave-number region. This resulted in absorption and scatteringcross-sections of 1.1× 10−7 cm2 and 5.9× 10−7 cm2, respectively (ω0

= 0.84). Figure 5.4 shows the single scattering albedo for spherical iceparticles in the mid-infrared as a function of wave-number and parti-cle radius. Obviously, the chosen values for ω0 cover a large fractionof the overall variability.

5.2 KOPRA - A single scattering model for the IR 99

1000 1500 2000 2500Wavenumber [cm-1]

0

10

20

30

40R

adiu

s [1

0-6 m

]

0.0

0.2

0.4

0.6

0.8

1.0

Figure 5.4: Single scattering albedo ω0 for spherical ice particles of differentradius covering the spectral region of the MIPAS/Envisat experiment. Forthe Mie calculations the refractive indexes of ice by Toon et al. Toon et al.(1994) were used.

Optical depths for the different cloud scenarios are given in Table 5.2for nadir direction and limb view with a tangent altitude at 11 km inthe middle of the cloud. In case of nadir geometry only scenario 5 isoptically thick, while in limb direction scenarios 3–5 are opaque.

Table 5.2: Cloud scenarios used for the ARTS/KOPRA intercomparison.Optical depths are given for the two cases ω0 = 0.24 and ω0 = 0.84 (inbrackets).

Cloud Number Volume Optical Opticalscenario density density depth depth

[cm−3] [µm3 cm−3] nadir limba

1 0.01 4.0194 1.68(1.75)×10−3 0.17(0.176)2 0.1 40.194 1.68(1.75)×10−2 1.7(1.76)3 1 401.94 1.68(1.75)×10−1 17.0(17.6)4 10 4019.4 1.68(1.75) 170(176)5 100 40194 16.8(17.5) 1700(1760)

a 11 km tangent altitude

100 5 Scattering model intercomparisons

5.2.3 Results for case ω0 = 0.24

Figure 5.5 shows the intercomparison between results from ARTS,and KOPRA(0)–KOPRA(3) for the case of strong absorption and fora line-of-sight crossing the middle of the cloud layer at 11 km tan-gent altitude. Spectra for lower altitudes are not shown here becausethe results were very similar. As reference, the top row compares thecloud-free model runs. In the region between the strong emission lineswhere the gaseous atmosphere is optically thin and to which, thus,the comparison between the cloud calculations will refer, differencesare below 0.5%.

The very thin cloud of scenario 1 mainly introduces a broadbandoffset above the clear sky spectrum due to the emission by the cloudparticles and the scattering of radiation from the troposphere and theearth’s surface. The scattered contribution is the difference betweenKOPRA(0) and ARTS and accounts for about 35% of the total ra-diance. Since the second zero scattering scheme KOPRA(1) modelsa larger emission from the cloud particles, the difference to ARTS iswith 15% somewhat reduced compared to KOPRA(0). The results ofboth single scattering models, KOPRA(2) and KOPRA(3) are nearlyidentical to the multiple scattering approach with less than 0.5% dif-ference.

In scenario 2, the radiance continuum is strongly increased andreaches with 2000 nW/(cm2 sr cm−1) the value of the Planck func-tion for the temperature at cloud altitude. The spectra of the singleand multiple scattering models clearly show signs of radiance of tro-pospheric origin, like the downward pointing absorption features ofthe water vapor lines. These structures are missing in the calcula-tions by both zero scattering schemes KOPRA(0) and KOPRA(1).The difference between those and ARTS are of the same magnitudeas in scenario 1. The accuracies of KOPRA(2) and KOPRA(3) withless then 1% are still much better in case of neglecting scattering.However, with less than 0.5% KOPRA(3) fits closer to ARTS thanKOPRA(2) with 1%.

In scenario 3 there is a further increase of the continuum radiancecompared to scenario 2. The background value for the scattering mod-

5.2 KOPRA - A single scattering model for the IR 101

0

1000

2000

3000

4000

H2O

CO2

CO2

H2O

CO2

CO2

CO2

Clear sky

0

1000

2000

3000

4000Scenario 1

0

1000

2000

3000

4000Scenario 2

0

1000

2000

3000

4000Scenario 3

0

1000

2000

3000

4000Scenario 4

946 947 948 949 950 951Wavenumber [cm-1]

01000

2000

3000

4000Scenario 5

KOPRA(3)KOPRA(2)KOPRA(1)KOPRA(0)ARTS

KOPRA(3)KOPRA(2)KOPRA(1)KOPRA(0)ARTS

Radi

ance

[nW

/(cm

2 sr cm

-1)]

-2

-1

0

-15

-10

-5

0

-15

-10

-5

0

-10

-5

0

-10

-5

0

5

946 947 948 949 950 951Wavenumber [cm-1]

-10-8-6-4-20

KOPRA(3) - ARTSKOPRA(2) - ARTSKOPRA(1) - ARTSKOPRA(0) - ARTS

KOPRA(3) - ARTSKOPRA(2) - ARTSKOPRA(1) - ARTSKOPRA(0) - ARTS

Diff

eren

ce [%

]

Figure 5.5: Comparison between limb spectra for 11 km tangent altitudecalculated with ARTS and different KOPRA options for the case of strongabsorption (ω0 = 0.24).

102 5 Scattering model intercomparisons

els is now larger than radiation of a blackbody at the cloud positionwhich is given by the result of KOPRA(1). The cloud is optically thickin limb, but not in nadir direction. Thus radiance from the warm tro-posphere and the earth’s surface increases the total signal. This isobvious from the fact that the spectra which include scattering, stillexhibit the downward pointing features of the tropospheric water ab-sorption lines. Differences show that the zero scattering models stillunderestimate the radiance by more than 10%. However, also the sin-gle scattering models differ from the ARTS reference by about 4.5%for KOPRA(2) and 2% for KOPRA(3).

The continuum radiance in scenario 4 is lower than in scenario 3due to the fact that, though the cloud is still not opaque in nadirdirection, less radiation from the troposphere reaches the particlesalong the line-of-sight which could scatter into the direction of the ob-server. Furthermore, the typical tropospheric absorption features arenot visible anymore. KOPRA(1) is with 3% difference closer to ARTSthan KOPRA(2) with about 9% difference. This is due to neglectingmultiple scattering. KOPRA(3) with 2.5% deviation from ARTS stilldelivers the best result.

In case of scenario 5, which is optically thick in nadir and limbdirection, the cloud closely resembles a black body. Thus, KOPRA(1)deviates from the multiple scattering model by only 0.1%. KOPRA(3)shows differences of 2% and KOPRA(2) of 11%.

5.2.4 Results for case ω0 = 0.84

Results of the model intercomparison for the case of large scattering(ω0 = 0.84) are shown in Figure 5.6. Though the optical depth of bothcases is not very different (see Table 5.2) the continuum signal in theoptically thin scenario 1 is increased by a factor of 1.4 compared tothe case of strong absorption. This is due to the increased scatteringof radiation originating in the warmer troposphere and on ground.Further, the flanks of the water vapor lines show downward pointingbroader features also caused by scattered tropospheric radiation. Withdifferences of 40% and more, the zero scattering models KOPRA(0)

5.2 KOPRA - A single scattering model for the IR 103

and KOPRA(1) are far off the reference while the single scatteringschemes KOPRA(2) and KOPRA(3) deviate by only 1.5% and 1%from ARTS.

With about 3000 nW/(cm2 sr cm−1) the background radiance of thereference spectrum for scenario 2 is 50% higher than the blackbodyradiance at cloud altitude. Very strong absorption features appearat the position of the water vapor and in the flanks of the CO2 lines.Thus, as in scenario 1, zero scattering models are not capable to modelthese effects. However, compared to scenario 2 of the strong absorptioncase, even single scattering models show large differences comparedto the ARTS reference: KOPRA(2) underestimates the radiance byabout 7% while KOPRA(3) calculates 3% lower values than ARTS.

For scenario 3, which is optically thick in limb direction, theradiance calculated by ARTS reaches its maximum at about3400 nW/(cm2sr cm−1). This is by a factor of 1.7 larger than theblackbody at cloud position. The single scattering model, however,does not follow this further increase, but results in lower radiancesthan in scenario 2. Thus, the differences with respect to the multiplescattering calculation increase up to 20% for KOPRA(3) and 35% forKOPRA(2).

For scenario 4, the ARTS radiances decreased and are only 1.2 timeshigher than the blackbody. Therefore, the zero scattering model KO-PRA(1) compares best, with differences of 13%. The single scatteringapproach of KOPRA(3) is far off with up to 32% lower radiances.

For scenario 5, which is optically thick in limb and nadir directions,ARTS and blackbody calculations (KOPRA(1)) are nearly identicalwith less than 1% difference. Thus, quasi no radiation from the lowertroposphere reaches the instrument any more. The comparison withKOPRA(3) is with 8% differences much better than for scenarios 3and 4.

5.2.5 Summary and discussion

The validity range of zero and single scattering calculations for simula-tions of mid-IR limb emission measurements of clouds was investigated

104 5 Scattering model intercomparisons

0

1000

2000

3000

4000

H2O

CO2

CO2

H2O

CO2

CO2

CO2

Clear sky

0

1000

2000

3000

4000Scenario 1

0

1000

2000

3000

4000Scenario 2

0

1000

2000

3000

4000Scenario 3

0

1000

2000

3000

4000Scenario 4

946 947 948 949 950 951Wavenumber [cm-1]

01000

2000

3000

4000Scenario 5

KOPRA(3)KOPRA(2)KOPRA(1)KOPRA(0)ARTS

KOPRA(3)KOPRA(2)KOPRA(1)KOPRA(0)ARTS

Radi

ance

[nW

/(cm

2 sr cm

-1)]

-2

-1

0

-40

-30

-20

-10

0

-40-30-20-10

010

-40-30-20-10

010

-80-60-40-20

020

946 947 948 949 950 951Wavenumber [cm-1]

-60

-40

-20

0

KOPRA(3) - ARTSKOPRA(2) - ARTSKOPRA(1) - ARTSKOPRA(0) - ARTS

KOPRA(3) - ARTSKOPRA(2) - ARTSKOPRA(1) - ARTSKOPRA(0) - ARTS

Diff

eren

ce [%

]

Figure 5.6: Comparison between limb spectra for 11 km tangent altitudecalculated with ARTS and different KOPRA options for the case of strongscattering (ω0 = 0.84).

5.2 KOPRA - A single scattering model for the IR 105

by comparison with the multiple scattering model ARTS-DOIT. Sce-narios from optically thin to thick clouds with a low (ω0 = 0.24) anda high (ω0 = 0.84) single scattering albedo were used as baseline forthe calculations.

For cloud scenarios 1 and 2, which are optically thin in limb di-rection, the single scattering approaches achieve results with maxi-mum errors of a few percent. These cloud scenarios resemble polarstratospheric clouds and subvisible cirrus clouds. Thus single scatter-ing models are sufficient for evaluation of such measurements. Zeroscattering schemes, however, show errors of more than 15% and 40%,depending on the single scattering albedo. From Figure 5.4 it is clearthat the zero scattering approaches can only be used for particle radiiless then about 1 µm at the lower and less than 0.2 µm at the higherend of the shown spectral interval.

For clouds, which are optically thick in limb direction, the qualityof single scattering calculations strongly depends on the single scat-tering albedo. The differences for the model KOPRA(3) range fromonly 2 – 3 % for ω0 = 0.24 up to 10 – 30% for ω0 = 0.84. For the lat-ter case larger errors appear for clouds which are optically thick inlimb, but not in nadir direction and which, thus, still scatter a largeamount of lower tropospheric radiation into the instrumental line ofsight. Combining these results with Figure 5.4 it is clear that for theregions with ω0 > 0.8, which are situated mainly above 1000 cm−1 forradii between 1 and 10–20 µm, the single scattering approaches deliverresults with uncertainties in the range of the investigated high scat-tering case. However, in the atmospheric window below 1000 cm−1,where the ice absorption peak is located, single scattering simulationsin case of optically thick clouds should be reliable for particle sizes ofless than 10–20 µm. For larger particles, the single scattering albedo isaround 0.5 over the whole wavelength region. Here, we estimate thatthe accuracy of single scattering calculations lies between the extremecases and is in the range of about 5 – 15%.

Comparing the results of different KOPRA options, we see, thatKOPRA(1) and KORA(3) achieve better results than KOPRA(0) andKOPRA(2). The latter implementations underestimate radiances be-cause they neglect the scattering source term while KOPRA(1) and

106 5 Scattering model intercomparisons

KOPRA(3) compensate for this by increasing the locally emitted ra-diation by replacing the absorption coefficient by the extinction coef-ficient in the direct and the scattered rays, respectively.

A further outcome of this study is that in cases with large singlescattering albedo and clouds which are optically thick in limb, butthin in nadir direction, the continuum radiance can be by a factor ofup to 1.7 enhanced with respect to the blackbody radiation at cloudtop altitude. This is a consequence of the scattering of radiation fromthe warm troposphere and the earth’s surface into the line-of-sight ofa limb viewing instrument. Such strong effects are expected to be de-tectable in recently measured spectra by MIPAS on Envisat or in dataof previous missions of mid-IR limb emission sounders like CRISTA(Cryogenic Spectrometers and Telescopes for the Atmosphere) (Spanget al., 2001) or CLEAS (Cryogenic Limb Array Etalon Spectrometer)(Roche et al., 1993).

5.3 Comparison between the Monte Carloand the Discrete Ordinate approach

This section describes the comparison between the two ARTS inter-nal scattering algorithms. In contrast to the previous sections, wherethe 1D DOIT algorithm was compared to other 1D scalar models,this section presents a comparison of 3D fully polarized models. TheMonte Carlo and the DOIT models are the first models of this kindfor radiative transfer modeling in the microwave region.

5.3.1 The Monte Carlo approach

A reversed Monte Carlo method has been implemented as a secondscattering module besides the DOIT module into the ARTS modelby Cory Davis. A strong consideration here was that the simplicityof the Monte Carlo radiative transfer concept should translate to re-duced development time. Also, reversed Monte Carlo methods allowall computational effort to be concentrated on calculating radiances

5.3 ARTS - Monte Carlo Scattering Model 107

for the desired line of sight, and the nature of Monte Carlo algorithmsmakes parallel computing trivial.

Among the available Backward Monte Carlo RT models, severaldo not allow a thermal source, or do not consider polarization fully,which means that they can not handle a non-diagonal extinction ma-trix (e.g., Liu et al. (1996)). Some consider neither thermal sourcenor polarization (e.g., Oikarinen et al. (1999), Ishimoto and Masuda(2002)).

The ARTS Monte Carlo algorithm is described in Davis et al.(2004). The flowchart shown in Figure 5.7 illustrates the algorithm.In the following a short summary will be given without looking intothe details:1. Begin with a new photon at the cloud box exit point and sample

a path length ∆s along the first requested line of sight using theprobability density function g0(∆s), which depends on the evolutionoperator O and on k, which is related to the extinction matrix.

2. A random number r is drawn to choose between emission and scat-tering using a quantity similar to the scattering albedo, ω.

3. If ω < r, the event is considered to be emission, the reversed raytracing is terminated and the Stokes vector contribution Ii is calcu-lated using the Planck function Ib(T ), where T is the temperature.Then return to step 1.

4. If ω > r, the event is considered to be scattering. At the scatteringpoint a new incident direction (θinc, φinc) is sampled according tothe probability function g(θinc, φinc), which is calculated using thephase matrix and the extinction matrix. The contribution to theintensity from this incident direction is obtained by the operatorQk.

5. A path length is sampled along the new direction.6. If this path length leads the photon outside the cloud box the con-

tribution of this photon on the cloud box boundary Ii is calculatedusing the matrix Qk.

7. Otherwise, if the sampled path length keeps the path within thecloud box, return to step 3.

8. When the Nth photon has reached the boundary of the cloud boxthe Monte Carlo algorithm is stopped and the field on the boundary

108 5 Scattering model intercomparisons

of the cloud box is used as the radiative background for the finalcloud-sensor clear sky radiative transfer.

5.3.2 Setup

Atmosphere

The comparison simulations were performed for two frequencies, 122and 230 GHz. These frequencies correspond to channels of the EOSMLS instrument, which can be used for cloud studies. Simulations forthis instrument are shown in Chapter 8. Atmospheric profiles weretaken from the FASCOD (Anderson et al., 1986) data for tropicalregions. For 122 GHz, only the species O2, N2 and H2O needed tobe included. For 230 GHz, CO and O3 were added. This selectionof gaseous species is based on Waters et al. (1999). The water vaporprofile was adjusted so that the relative humidity is 100% with respectto ice at altitudes with non-zero ice mass content.

Cloud scenario

A thin cirrus cloud layer was selected. The most simple case to startwith is a box-shaped cloud in pressure, latitude and longitude coor-dinates. Since the Monte Carlo method is implemented only for 3Dclouds, a 1D comparison is not possible. The FASCOD profiles how-ever are 1D, so a box-shaped 3D cloud was embedded into the 1D at-mosphere. The particle size distribution by Mc Farquhar and Heyms-field, which has been described in Section 3.5.3, was used. A singleaspect ratio of 1.5 was applied for the whole cloud. It was assumedthat the cloud particles are horizontally aligned. The ice mass contentof the cloud was 0.1 g/m3 and the cloud altitude was 11.9 – 13.4 km.The horizontal extent was 400 km×400 km, which corresponds to 3.6

latitude times 3.6 longitude for tropical regions. Two different sets oflines of sights were considered. They are illustrated in Figure 5.8. Alllines of sight of set A intersect in the cloud at a point located 50 kmaway from the north edge. The lines of sight of set B intersect in thecloud in the center. Clear sky and cloudy simulations were performed

5.3 ARTS - Monte Carlo Scattering Model 109

SCAT

TERIN

G

sam

ple

anew

incid

ent

directio

n(θ

inc ,φ

inc )

ac-

cord

ing

to

g(θ

inc ,φ

inc )

=Z

11 (θ

scat ,φ

scat ,θ

inc ,φ

inc )

sin(θ

inc )

K11 (θ

scat ,φ

scat )−

Ka1 (θ

scat ,φ

scat )

Calcu

late

the

matrix

Qk

=Q

k−

1 qk

,w

here

qk

=sin

(θinc )

k O(s

k,s

k−

1 )Z(n

k−

1 ,nk )

g(∆

s)g(θ

inc ,φ

inc )ω

,

and

Q0

=

.

Beg

inat

the

cloud

box

exitpoin

tw

itha

new

photo

n.

Sam

ple

apath

length

,∆

salo

ng

the

first

line

ofsig

ht

usin

gth

ePD

F

g0 (∆

s)=

kO

11 (∆

s)

1−

O11 (u

0 ,s0 )

.

START

i=

N?

CLO

UD

EXIT

ST

OK

ES

VECT

OR

I(n,s

0 )=

O(u

0 ,s0 )I(n

,u0 )

+〈I

i(n,s

0 )〉.

Use

this

as

the

radia

tiveback

gro

und

for

final

cloud-sen

sor

clearsk

yRT

.

FIN

ISH

EM

ISSIO

N

Ii(n

,s0 )

=Q

kO

(sk+

1 ,sk )K

a (nk,s

k+

1 )Ib (T

,sk+

1 )

g(∆

s)(1

−ω

)

BO

UN

DARY

Ii(n

,s0 )

=Q

kO

(uk,s

k )I(nk,u

k )

O11 (u

k,s

k )

EM

ISSIO

N

Ii(n

,s0 )

=O

(s1 ,s

0 )Ka (n

0 ,s1 )I

b (T,s

i )

g0 (∆

s)(1

−ω

)

Sam

ple

anew

path

length

,∆

salo

ng

the

new

directio

nusin

gth

ePD

F

g(∆

s)=

kO

11 (∆

s)

OU

TSID

E?

r>

ω?

r>

ω?

NO

NO

NO

YES

YES

YES

YES

i=

1

k=

0

NO

i=

i+

1

k=

k+

1

Figure 5.7: Flowchart illustrating the Monte Carlo algorithm. Courtesy ofCory Davis.

110 5 Scattering model intercomparisons

for tangent altitudes between 1 and 13 km. Since the DOIT modelyields the full radiation field at the same time, a fine tangent altituderesolution was used (100 m). Monte Carlo simulations were performedon a tangent altitude grid with a resolution of 1 km. The zenith an-gle grid for the DOIT calculation was optimized to a grid accuracyof 0.1%. In 3D, so far only linear interpolation is implemented. Thismeans that at critical zenith angles the error in the calculation can goup to approximately 1% (compare Figure 4.7). For most altitudes itis expected to be less than 0.2%. The error depends significantly onthe cloud optical thickness. For very thick clouds it can be larger thanthe estimated result, since more iterations are required.

5.3.3 Results

The results for 122 GHz are presented in Figure 5.9. The top panelshows the obtained radiances in Rayleigh Jeans BT for tangent alti-tudes from 1 to 13 km. The grey line is the clear sky field. Black linescorrespond to DOIT results and the markers are the Monte Carloresults for both LOS sets. In the radiance plot, there are no obviousdifferences between the DOIT and the Monte Carlo results or betweenthe LOS set A and LOS set B. The middle left panel shows the dif-ference between cloudy and clear sky radiances. Here we see that theMonte Carlo model shows less BT depression than the DOIT model.The absolute difference, which is shown in the bottom left panel, isbetween 0.2 and 0.4 K. The middle right panel shows the results ob-tained for the polarization difference Q, also in Rayleigh Jeans BT.Here also the Monte Carlo model shows slightly smaller values. Theabsolute difference between the models is approximately 0.05 K. Forboth, I and Q, the sign of the differences between the models is at12 km tangent altitude opposite to the sign of the differences at othertangent altitudes. This can be explained by two opposing mechanisms:scattering away from the LOS for low tangent altitudes and scatteringinto the LOS for high tangent altitudes. If the scattering effect in theDOIT model was slightly greater than in the Monte Carlo model, theresult would show more BT depression for low tangent altitudes and

5.3 ARTS - Monte Carlo Scattering Model 111

A

B

−1.8

0.0

1.8

lat.

viewing direction

−1.8 0.0 1.8lon.

-200 -100 0 100 200

horizontal extent [km]

6376

6378

6380

6382

6384

pre

ssu

re [

hP

a]

A

B

A

B

Figure 5.8: Top: Cloud box and crossing points of LOS sets A and B asseen from the top. The arrow shows the viewing direction of the instrument.Bottom: Lines of sight for tangent altitudes 1 – 13 km of sets A and B andthe cloud box.

112 5 Scattering model intercomparisons

more BT enhancement for high tangent altitudes. At 12 km tangentaltitude, scattering away from the LOS still dominates, but the BTdepression becomes smaller. The middle left panel shows irregularitiesin the DOIT result due to interpolation. These would disappear witha finer zenith angle grid. The differences between the DOIT and theMonte Carlo models are similar for both sets of LOS, A and B. Theagreement between the models is very good. It is well inside the errorestimates of the DOIT model.

The results for 230 GHz are presented in Figure 5.10. The panelsare arranged in the same way as in Figure 5.9. The radiance plotshows that the scattering effect is much larger for 230 GHz comparedto 122 GHz. There is a BT depression of approximately 50 K for tan-gent altitudes below 8 km and a BT enhancement of up to 100 K fortangent altitudes above 8 km. At approximately 8 km there is almostno difference between cloudy and clear sky radiances. That means thatat this point the same amount of radiation is scattered into the LOSas away from the LOS. The polarization difference is also much largercompared to 122 GHz, it goes up to 8 K. The middle panels show noobvious deviations between the Monte Carlo model and the DOITmodel, but the bottom panels show, that the difference is almost 2 Kfor small tangent heights and 4 K at 12 km tangent altitude. For Q

the absolute differences are small apart from 12 km tangent altitudewhere it goes up to -2 K. One of the reasons for this deviation is theinterpolation of the radiation field. Since many iterations are requiredfor the rather strong cloud the interpolation error could be more than1%.

5.3.4 Discussion

At 122 GHz the Monte Carlo model and the 3D DOIT model are invery good agreement, the differences are less than 0.5 K for total inten-sities and about 0.05 K for the polarization difference. At 230 GHz thedifference is much larger, for most tangent altitudes approximately 2 Kfor the intensities and approximately 0.3 K for the polarization differ-ence. The major differences between the two frequencies are the gas

5.3 ARTS - Monte Carlo Scattering Model 113

Monte Carlo (MC)vs. DOIT

Frequency:122 GHz

190 200 210 220 2300

5

10

15

I [ K ]

z tan [

km ]

A DOITA MCB DOITB MCclear sky

−4 −3 −2 −1 00

5

10

15

Icloudy

− Iclear

[ K ]

z tan [

km ]

0 0.2 0.4 0.6 0.80

5

10

15

Q [ K ]

z tan [

km ]

−0.4 −0.2 0 0.2 0.40

5

10

15

∆ I [ K ]

z tan [

km ]

−0.1 −0.05 0 0.05 0.10

5

10

15

∆ Q [ K ]

z tan [

km ]

A DOITA MCB DOITB MC

A DOIT−MCB DOIT−MC

Figure 5.9: Results obtained for 122 GHz: The top panel shows the abso-lute radiances and the clear sky radiances (grey line). The middle left panelshows the intensity differences and the middle right panel shows the po-larization signal. The bottom panels show the absolute differences betweenthe DOIT and the MC modules for intensity and polarization.

114 5 Scattering model intercomparisons

Monte Carlo (MC)vs. DOIT

Frequency:230 GHz

0 100 200 3000

5

10

15

I [ K ]

z tan [

km ]

A DOITA MCB DOITB MCclear sky

−100 −50 0 50 1000

5

10

15

Icloudy

− Iclear

[ K ]

z tan [

km ]

−10 −5 0 5 100

5

10

15

Q [ K ]

z tan [

km ]

−6 −4 −2 0 20

5

10

15

∆ I [ K ]

z tan [

km ]

−2 −1 0 10

5

10

15

∆ Q [ K ]

z tan [

km ]

A DOITA MCB DOITB MC

A DOIT−MCB DOIT−MC

Figure 5.10: Results obtained for 230 GHz: The top panel shows the abso-lute radiances and the clear sky radiances (grey line). The middle left panelshows the intensity differences and the middle right panel shows the po-larization signal. The bottom panels show the absolute differences betweenthe DOIT and the Monte Carlo modules for intensity and polarization.

5.3 ARTS - Monte Carlo Scattering Model 115

absorption and the single scattering properties. The clear sky fieldshows, that the gas absorption at high altitudes is much larger at122 GHz compared to 230 GHz. The single scattering properties in-crease with frequency as can be seen from Figure 3.6. Therefore, cloudscattering has a much higher impact at 230 GHz.

The comparison shows, that both models are able to simulate theeffect of cirrus clouds. The features of the scattering signal, i.e., onlyBT depression at 122 GHz and BT depression for low tangent alti-tudes and BT enhancement for high tangent altitudes at 230 GHz, arewell brought out in the results. The numerical approaches to solve theVRTE Equation (1.42) are completely different, hence different nu-merical errors are involved. The accuracy of the Monte Carlo methodis mainly restricted by the number of photons used for the calcula-tions. In principle an arbitrary accuracy can be obtained by increasingthe number of photons. The drawback is an increase in computationtime.

A disadvantage of the 3D DOIT model is, that the numerical gridscan not be chosen arbitrarily fine. The required memory for the com-putation depends on the size of the radiation field, which is dis-cretized in pressure, latitude, longitude, incoming and scattered di-rections (compare Section 4.1.1). Also the computation time dependsmainly on the discretization of the radiation field. Since the zenithangle grid is used for incoming and scattered directions, a doublingof the number of zenith angle grid points leads to four times greatercomputation memory and time requirements. A very fine zenith anglegrid is unavoidable when one wants to achieve accurate results, sincethe radiation field strongly increases at zenith angels of approximately90. Doubling the size of the cloud box in all three spatial coordinateseven leads to eight times greater computation time and memory re-quirements. The zenith angle grid discretization is not essential in theMonte Carlo model, since the radiation field is not interpolated in thisdimension. Computation time also increases with the size of the cloudbox, but it does not factorize like in the DOIT model. Another nu-merical error source is the calculation of 3D propagation paths. Prop-agation paths in a 3D atmosphere are calculated using an iterativeapproach, in contrast to 1D, where they are calculated analytically.

116 5 Scattering model intercomparisons

When the number of iterations increases, the numerical errors due tothese calculations also increase.

Qualitatively the differences between the DOIT model and theMonte Carlo model can be understood, but the comparison showsthat it is difficult to make correct error estimates. In both models,for 230 GHz, the estimated accuracy was better than the differencebetween the models, which means that there are additional numericalerrors which need to be taken into account.

A feature of the DOIT model is that it yields the whole radiationfield. This can be useful to obtain a basic physical understanding of forexample polarization effects in the cloud. However, for limb sounding,when one is only interested in a few selected viewing directions, theMonte Carlo model should be preferred presently, if the cloud fieldis strongly inhomogeneous, so that a fine spatial discretization is re-quired. In the future, with faster processors and more computationmemory, one could also use the DOIT module for bigger cloud sce-narios. For a thin cirrus layer with a large horizontal extent, it makessense to use the 1D DOIT model, which is much faster than the MonteCarlo model.

5.3.5 Summary and conclusions

This comparison study was very important, since the two comparedmodels are the very first ones, which can be used to simulate limbradiances in 3D spherical atmospheres for the microwave region. Atthe moment no other models are available for comparison. The resultof this study is very promising. The agreement between the modelsis satisfactory. It shows that both the Monte Carlo method and thediscrete ordinate method can be applied for solving the VRTE in a 3Dspherical atmosphere. Optimization in accuracy and speed is plannedto be implemented in the 3D-DOIT model in the future.

6 Unpolarized 1D simulations tostudy the impact of cirrus cloudson microwave limb measurementsof the MASTER instrument

In this chapter the first simulations of microwave limb radiances withclouds are presented and analysed. They are computed using the 1Dunpolarized version of the DOIT scattering module. Limb spectraare generated for the frequency bands of the MASTER (MillimeterWave Acquisitions for Stratosphere/Troposphere Exchange Research)instrument (Buehler, 1999). The impact of various cloud parametersis investigated. Simulated brightness temperatures most strongly de-pend on particle size, ice mass content and cloud altitude. The impactof particle shape is much smaller, but still significant. Increasing theice mass content has a similar effect as increasing the particle size,this complicates the prediction of the impact of clouds on microwaveradiances without exact knowledge of these parameters. The workpresented in this chapter has been published in Emde et al. (2004b).

6.1 General setup for the simulationsBefore defining the cloud parameters the general setup of the modelneeds to be defined. This includes the composition and the geometryof the model atmosphere, the sensor and the numerical setup.

117

118 6 Microwave limb spectra

200 250 3000

10

20

30

40

50

60T

T [K]

z [k

m]

10−8

10−6

0

10

20

30

40

50

60O3

VMR [ ]

z [k

m]

10−6

10−5

10−4

10−3

10−2

10−1

0

10

20

30

40

50

60H2O

VMR [ ]

z [k

m]

Figure 6.1: Atmospheric profiles for temperature, ozone and water vaportaken from FASCOD data for mid-latitudes in summer.

6.1.1 Atmosphere

The simulations are performed in a 1D spherical model atmosphere.Of course a 1D atmosphere is not a realistic environment for model-ing clouds, because clouds are horizontally strongly inhomogeneous.Each cloud included in a 1D model corresponds to full cloud coveragearound the globe. Thus the 1D calculations presented in this study canonly be taken as an upper limit of scattering effects on the simulatedradiances. Besides nitrogen and oxygen the two major atmosphericgases, water vapor and ozone, are included. The concentrations aretaken from FASCOD (Anderson et al., 1986) data for mid-latitudes insummer. Profiles for temperature, ozone and water vapor are shownin Figure 6.1. Gas absorption is calculated based on the HITRAN(Rothman et al., 1998) molecular spectroscopic database using theARTS model (first version ARTS-1-0). Refraction has been neglectedin all calculations.

6.1 General setup for the simulations 119

6.1.2 Sensor setup

The spectral ranges of bands B (293 – 306 GHz), C (317 – 326 GHz),D (342 – 349 GHz) and E (496 – 506 GHz) of the MASTER instrumentare used for the simulations. The different scans correspond to tangentaltitudes from 0 km up to 12.5 km, where 0 km altitude is the Earth’ssurface and 12.5 km is a tangent altitude above the cloud. For simplic-ity, the sensor is assumed to be ideal, which means that it measuresexactly the intensity of the incoming radiation and is not subject tonoise or other errors. This implies that the signal neither depends onthe antenna pattern nor on the polarization of the measured radia-tion. The scalar version of the DOIT module is applied, thus only thefirst component of the Stokes vector is calculated, this means thatpolarization is neglected in the radiative transfer.

6.1.3 Numerical setup

For the accuracy of the results the discretization, especially in the an-gular domain, is very important. For the calculation of the scatteringintegral Equation (4.5), an equidistant zenith angle grid with an incre-ment of 10 was taken. This is not sufficient for the radiative transfercalculation (Equation (4.7)), since the radiation field is strongly inho-mogeneous around 90. This problem and its solution are describedin Section 4.3. For this study the zenith angle grid was optimized torepresent the clear sky radiation field with an accuracy of 0.1%. Thevertical grid is equidistant in altitude and the grid step-size is 0.5 km.If a propagation path step, i.e., the distance of two successive intersec-tion points of the line of sight (LOS) with the vertical grid, is longerthan 1 km, which occurs in limb geometry close to the tangent point,this step is divided into smaller (< 1 km) equidistant steps. The errorof the results is estimated to be in the range of 0.5% including theinterpolation of the radiation field as the main source for numericalerrors.

120 6 Microwave limb spectra

6.2 Definition of cloud scenarios for theinvestigation of the impact of differentcloud properties on limb radiances

To study the impact of different cloud properties (IMC, particle size,particle shape, altitude and frequency), the test cases compiled inTable 6.1 were defined. In all test cases, except case (3), it was assumedthat the cloud consists only of spherical particles. Band C is taken forall calculations except for case (5). As the scattering signal dependsalso on the amount of gas absorption, it is interesting to study theeffects in a band where we find frequency regions with high and otherswith low gas absorption (window regions). In this sense band C is thebest selection.

Cloud height 10 – 12 km means that the scattering properties are de-fined on all pressure grid points between 10 and 12 km. The propertiesare linearly interpolated between the grid points, in this case they areinterpolated between 9.5 and 10 km and between 12 and 12.5 km. Be-tween 10 and 12 km the scattering properties are constant. the gammadistribution, which was described in Section 3.5.2, was used to calcu-late the single scattering properties and the particle number densityfields.

Table 6.1 includes the following scenarios:

1. To investigate the impact of particle size spectra for the MASTER-C frequency band were calculated. The IMC is constant in all cal-culations (1.6·10−3 g/m3) and the mean effective radius varies from21.5 to 128.5 µm.

2. The influence of cloud altitude is also studied using band C. IMCand effective radius are constant, 1.6·10−3 g/m3 and 34.0 µm respec-tively. Calculations are performed for three different cloud altitudes:6 – 8 km, 8 – 10 km and 10 – 12 km.

3. This is the only case were the cloud is assumed to consist of cylin-drical particles. Spectra for five different aspect ratios in the rangefrom 0.3 to 4.0 are calculated. The aspect ratio is the diameterof the cylinder divided by its length. Again band C is used and

6.3 Results 121

IMC and effective radius are constant, 1.6·10−3 g/m3 and 68.5 µmrespectively.

4. Measurements have shown that there is a correlation between par-ticle size and ice mass content. Realistic scenarios according to aplot shown in Evans et al. (1998), which includes results from theFIRE (Kinne et al., 1997) campaign, were picked out. The IMC isvaried from 4·10−5 to 0.04 g/m3 and the effective radius from 21.5to 128.5 µm. Like in the other cases band C is taken and the cloudaltitude is 10 – 12 km.

5. Calculations for different frequency bands, B, C, D and E are per-formed to see the impact of the same cloud in different frequencyregions. IMC and effective radius are constant, 1.6·10−3 g/m3 and34.0 µm respectively.

6.3 Results6.3.1 Impact of particle size

Simulations for varying particle sizes and constant IMC are shown inFigure 6.2. The simulated radiances are presented in Rayleigh-Jeansbrightness temperature (BT) units. On the left hand side we see theresults at 8 km tangent altitude. The spectrum of the cloud consistingof the smallest particles (Reff = 21.5 µm) is almost identical to theclear sky spectrum. All other clouds show a BT depression in the win-dow regions. For a particle size of 34.0 µm the brightness temperaturedepression (∆BT) is about 4 K. It becomes much larger with increas-ing particle size. For a particle size of 68.5 µm, ∆BT is approximately17 K, for 85.5 µm approximately 28 K and for 128.5 µm approximately52 K.

The plots on the right hand side show the results at 11.5 km tangentaltitude, which look very different from those at 8 km tangent altitude.In the window regions about 317 and 322.5 GHz a BT enhancement isobserved. This means, that more radiation is scattered into the line ofsight (LOS) than away from the LOS. The largest BT enhancementcan be observed at about 318 GHz in the window region of band C.

122 6 Microwave limb spectra

Table 6.1: Definition of five test cases to study the effect of cloud propertieson limb radiances

IMC[g/m3]

Reff

[µm]cloud alt.[km]

band aspectratio

1 1.6·10−3 21.534.068.585.5128.5

10 – 12 C —

2 1.6·10−3 34.0 6 – 88 – 1010 – 12

C —

3 1.6·10−3 68.5 10 – 12 C 0.3, 0.5,1.0, 2.0,4.0

4 4·10−5 21.5 10 – 12 C —1.6·10−3 34.08·10−3 68.50.016 85.50.04 128.5

5 1.6·10−3 34.0 10 – 12 B, C,D, E

—

It ranges from about 5 K for a particle size of 21.5 µm and goes upto about 50 K for a particle size of 128.5 µm. Above approximately322.5 GHz we see a BT depression for all particle sizes. The totalextinction coefficient consists of particle absorption, particle scatteringaway from the LOS and gas absorption. If the contribution from thegas absorption is dominant, the scattering effect becomes small. Moreradiation is absorbed than scattered into the LOS. The BT depressionvaries from almost 0 K for the cloud with the smallest particles toapproximately 35 K for the largest particles. The results show, thatthe size of the particles has a very large impact on limb radiances.

6.3 Results 123

318 320 322 324

180

190

200

210

220

230

240

Frequency [ GHz ]

BT

[ K ]

8 km tangent height

318 320 322 324

−50

−40

−30

−20

−10

Frequency [ GHz ]

∆ B

T [ K

]

318 320 322 324100

120

140

160

180

200

220

240

Frequency [ GHz ]

BT

[ K ]

11.5 km tangent height

318 320 322 324

−30

−20

−10

0

10

20

30

40

Frequency [ GHz ]

∆ TB

[ K

]

Figure 6.2: Impact of particle size. Left panel – limb spectrum at 8 kmtangent altitude. Right panel – limb spectrum at 11.5 km tangent altitude.Clear sky spectrum (grey), and cloudy spectra for Reff = 21.5 µm ( — ),Reff = 34.0 µm (– · –), Reff = 68.5 µm (· · ·), Reff = 85.5 µm (– – –) andReff = 128.5 µm ( — ). The top plots show absolute BTs, the bottom plotsshow the differences between cloudy and clear sky spectra.

6.3.2 Impact of cloud altitude

The impact of cloud altitude on limb radiances is presented in Fig-ure 6.3. The IMC was the same as taken for Figure 6.2. The assumedparticle size for those calculations was 34.0 µm. This is one of the op-tically thinner clouds. The grey line is the clear sky spectrum. Wesee that the cloud at altitude 6 – 8 km has a very small impact on theradiances, below 0.5 K at 318 GHz and even less at higher frequencies.The cloud at 8 – 10 km altitude leads to a BT depression of maximal1.7 K and the one at 10 – 12 km to a BT depression of maximal 2.4 K.

On the right side the BT is plotted as a function of tangent altitude.The lowest cloud shows in all tangent altitudes only a very smallimpact. The cloud at 8 – 10 km altitude leads to a BT depression attangent altitudes up to 9 km and then to small BT enhancement up to

124 6 Microwave limb spectra

318 320 322 324

215

220

225

230

235

240

Frequency [ GHz ]

BT

[ K ]

7 km tangent height

318 320 322 324

−2

−1.5

−1

−0.5

0

Frequency [ GHz ]

∆ B

T [ K

]

−2 0 2 4 6 8 100

2

4

6

8

10

12

BT [ K ]

tang

ent h

eigt

[ km

]

318 GHz

Figure 6.3: Impact of cloud altitude. Left panel - limb spectrum at 7 kmtangent altitude. Clear sky spectrum (grey), and cloudy spectra for cloudaltitude 6 – 8 km ( — ), 8 – 10 km (– · –) and 10 – 12 km (· · ·). Top plotshows absolute BTs, bottom plot differences from the clear sky case. Rightpanel - radiances at 318 GHz as a function of cloud altitude. 6 – 8 km ( — ),8 – 10 km (– · –) and 10 – 12 km (· · ·).

10.5 km tangent altitude. The cloud signal is observed up to 10.5 km,since the scattering properties are interpolated linearly between thegrid points. As mentioned above, the cloud ranges from 7.5 to 10.5 kmin this case. For the highest cloud, we also observe a BT depression upto 9 km tangent altitude and a very high BT enhancement from 10 to12 km tangent altitude. The enormously higher BT enhancement forthe high cloud is due to the fact that in lower altitudes water vaporabsorption is very large.

6.3.3 Impact of particle shape

The results of limb spectra for clouds consisting of cylindrical particleswith different aspect ratios varying from 0.3 to 4.0 is shown in Fig-ure 6.4. The calculations were done for the same IMC as taken in the

6.3 Results 125

318 320 322 324

210

215

220

225

230

235

240

Frequency [ GHz ]

BT

[ K ]

8 km tangent height

318 320 322 324

−15

−10

−5

Frequency [ GHz ]

∆ B

T [ K

]

318 320 322 324100

120

140

160

180

200

220

240

Frequency [ GHz ]

BT

[ K ]

11.5 km tangent height

318 320 322 324

−10

0

10

20

30

Frequency [ GHz ]

∆ B

T [ K

]

Figure 6.4: Impact of particle shape. Left panel - limb spectrum at 8 kmtangent altitude. Right panel - limb spectrum at 11.5 km tangent altitude.Clear sky spectrum (grey), and cloudy spectra for aspect ratios 0.3 ( — ),0.5 (– · –), 1.0 (· · ·), 2.0 (– – –) and 4.0 ( — ). The top plots show absoluteBTs, the bottom plots show the differences between cloudy and clear skyspectra.

previous calculations, but here a larger mean particle size of 68.5 µmwas taken. We can see immediately, that the difference between thecurves is small compared to the scattering effect itself. At 8 km tangentaltitude the total BT depression is about 17 K. The difference betweenthe spectra for different aspect ratios is only approximately 2 K. At11.5 km tangent altitude the maximal BT enhancement is about 30 Kand the differences for different shapes are again in the order of 2 K.At this point it does not make sense to interpret the results in detail,for example, whether plates (aspect ratio > 1) show a higher signalthan cylinders (aspect ratio < 1). It may only be concluded, that theeffect of particle shape is much smaller than the effects of particle sizeand cloud altitude.

126 6 Microwave limb spectra

6.3.4 Cloud scenarios with correlation betweenIMC and Reff

Results of the calculations for the “realistic” clouds are shown inFigure 6.5. The left side shows the results for 8 km tangent alti-tude. The impact of the optical thinnest cloud (IMC = 4·10−5 g/m3,Reff = 21.5 µm) is very small, the spectrum can not be distinguishedfrom the clear sky spectrum. The impact of the second cloud (IMC= 1.6·10−3 g/m3, Reff = 34.0 µm), which is the one used for studyingthe effect of cloud altitude is very small compared to the other cloudscenarios. For the optically thickest cloud the BT depression may goup to 120 K. The highest BT enhancement at 11.5 km is observed fora medium cloud thickness (IMC = 8·10−3 g/m3, Reff = 68.5 µm). Thereason is that in case of very thick clouds, mainly radiation scatteredinto the LOS in the upper part of the cloud contributes to the spectra.It is very probable that radiation scattered into the LOS in the loweror middle part of the cloud is again absorbed or scattered away fromthe LOS by a cloud particle. For thinner clouds, most radiation, whichis scattered once into the LOS, will continue to propagate into thisdirection without being disturbed by other cloud particles.

Although the micro-physical cloud properties, i.e., particle size andIMC are “realistic”, a limb instrument would observe smaller scat-tering signals. As mentioned above, these calculations can only betaken as an upper limit of scattering effect, because the 1D mode ofthe model was used, assuming a homogeneous cloud cover. This as-sumption is unrealistic, particularly for the intense cloud cases, whichcorrespond to clouds of limited horizontal extent.

6.3.5 Spectra for different frequency bands

In Figures 6.6 to 6.9 results for MASTER bands B, C, D and E arepresented. They correspond to cloud case 2. The spectra are shown forthe tangent altitudes of 3, 9, 10.5 and 11.5 km. The top panels showthe clear sky spectra, the middle panels the spectra in the presenceof clouds and the bottom panels show the difference ∆BT between

6.3 Results 127

318 320 322 324

120

140

160

180

200

220

240

Frequency [ GHz ]

BT

[ K ]

8 km tangent height

318 320 322 324

−100

−80

−60

−40

−20

Frequency [ GHz ]

∆ B

T [ K

]

318 320 322 324100

120

140

160

180

200

220

240

Frequency [ GHz ]

BT

[ K ]

11.5 km tangent height

318 320 322 324

−60

−40

−20

0

20

40

Frequency [ GHz ]

∆ B

T [ K

]

Figure 6.5: Impact of IMC. Left panel - limb spectrum at 8km tangent altitude. Right panel - limb spectrum at 11.5 km tan-gent altitude. Clear sky spectrum (grey), and cloudy spectra forIMC = 4·10−5 g/m3 ( — ), 1.6·10−3 g/m3 (– · –), 8·10−3 g/m3 (· · ·),0.016 g/m3 (– – –) and 0.04 g/m3 ( — ) and the corresponding particle sizes.The top plots show absolute BTs, the bottom plots show the differences be-tween cloudy and clear sky spectra.

cloudy and clear sky spectra. IMC, Reff and cloud altitude are definedin Table 6.1.

In band B there is an oxygen line at approximately 298.5 GHz. Theother visible spectral lines are due to ozone. For 3 km tangent altitudethere is a BT depression of maximal 1.5 K. All other tangent altitudesshow a BT enhancement due to clouds. The difference between cloudyand clear sky spectrum is very small (< 3 K) for a tangent altitude of9 km. This tangent altitude is below the cloud. As the considered cloudis rather optically thin, we can already see at 9 km a BT enhancement.For 10.5 and 11.5 km tangent altitude the BT enhancement can bemore than 15 K.

In band C there are ozone lines around 317 and 320 GHz. The line

128 6 Microwave limb spectra

at approximately 325 GHz is a water vapor line with a high absorp-tion. This means that the largest scattering effect can be observed inthe window regions around the ozone line at 320 GHz. Like in band Bwe see a BT enhancement for the tangent altitudes inside the cloud(10.5 km and 11.5 km). The maximum at 11.5 km tangent altitude fora frequency of about 318 GHz is with a brightness temperature differ-ence of 10 K smaller than the maximum in band B, because the totalgas absorption in band C is higher than in band B. The scatteringcoefficient increases with frequency in the microwave range (cf. Fig-ure 3.6). At 11.5 km tangent altitude (upper part of the cloud) mostradiation is scattered into the LOS. At 10.5 km tangent altitude, ∆BTis smaller compared to band B. As the scattering coefficient is largerin band C, the probability for multiple scattering is increased. If radi-ation is scattered into the LOS at 10.5 km there is a large probabilitythat it is again scattered out of the LOS during the propagation from10.5 to 12 km. The spectrum at tangent altitude 9 km shows alreadya BT depression in band C, the amount is less than 2 K. The BT de-pression at 3 km tangent altitude is maximal 2.5 K, about 1 K largerthan the maximal depression in band B. This is the expected result,as the extinction coefficient is increased for higher frequencies.

In band D there are ozone lines around 343.3 GHz and an oxygenline at approximately 345.3 GHz. The total gas absorption is similarto band B. Since the scattering coefficient is larger for higher frequen-cies, the BT enhancement from the cloud at 10.5 and 11.5 km tangentaltitude is increased, it is in this case maximal 11 K and 14 K respec-tively. At 9 km tangent altitude, the cloud effect is very small, onlya depression of less than 1 K is observed. The BT depression at 3 kmtangent altitude is approximately 3 K.

Band E is in a much higher frequency region with much largerscattering coefficients than the other bands. But also the total gasabsorption is the largest in this band. Only for a tangent altitude of11.5 km we see a BT enhancement throughout the whole band, themaximum value is about 17 K. The smallest brightness temperaturedifference in band E is observed for 10.5 km tangent altitude, whichcorresponds to an altitude inside the cloud. At 9 km tangent altitudethe BT depression is larger (max. 12 K) than at 3 km tangent altitude

6.3 Results 129

(max. 8 K). As the propagation path through the cloud is longer fora limb scan with 9 km tangent altitude than for a scan with 3 kmtangent altitude and the cloud extinction is large, more BT depressionis observed at a tangent altitude of 9 km.

6.3.6 Comparison with nadir radiances

Compared to nadir radiances, limb radiances are more complex. Innadir geometry, cirrus clouds always lead to a brightness temperaturedepression compared to the clear sky radiances. Nadir instrumentslook at the tropopause and the ground, where the major source ofmicrowave radiation is located. Clouds absorb and scatter part of theradiation out of the LOS. When we consider the lower atmosphere andthe emitting ground as major sources of radiation, the propagationdirection of the radiation is upwards. In nadir geometry we measurethe up-welling radiation, there can not be an enhancement in thisdirection due to scattering according to the law of energy conservation.

In limb geometry, clouds usually lead to a brightness temperaturedepression if the tangent altitude lies below the cloud. But is can alsolead to a BT enhancement if the tangent altitude is inside the cloud,because part of the up-welling radiation from the Earth’s surface andthe lower atmosphere is scattered into the LOS. In the clear sky casethe sensor does not see thermal emission from the lower atmosphereat all.

The impact of cloud in nadir geometry is smaller, as the path-lengththrough the cloud is much shorter. Of course the path-length throughthe cloud for a limb measurement depends on the horizontal extentof the cloud. Since clouds have an infinite extent in a 1D model, theradiances are overestimated for most cases. In reality the cloud cover-age is horizontally inhomogeneous. But still the path-length throughthe clouds is larger in limb. The clouds presented in this study aremostly optically thin. The highest brightness temperature depressionin nadir observing geometry by the clouds used to study the impact ofparticle size, case (1), is only 1.7 K. Only case (4), where particle size

130 6 Microwave limb spectra

294 296 298 300 302 3040

50

100

150

200

250

Frequency [ GHz ]

BT

[ K ]

Clear

h=3 kmh=9 kmh=10.5 kmh=11.5 km

294 296 298 300 302 3040

50

100

150

200

250

Frequency [ GHz ]

BT

[ K ]

Cloudy

h=3 kmh=9 kmh=10.5 kmh=11.5 km

294 296 298 300 302 304

0

5

10

15

Frequency [ GHz ]

∆ B

T [ K

]

Cloudy − Clear

h=3 kmh=9 kmh=10.5 kmh=11.5 km

Figure 6.6: Limb spectra at different tangent altitudes for MASTER bandB (294 – 304 GHz). The top panel shows the clear sky spectra, the middlepanel shows the cloudy spectra and the bottom panel shows the differencebetween cloudy and clear sky spectra.

6.3 Results 131

317 318 319 320 321 322 323 324 325

100

150

200

250

Frequency [ GHz ]

BT

[ K ]

Clear

h=3 kmh=9 kmh=10.5 kmh=11.5 km

317 318 319 320 321 322 323 324 325

100

150

200

250

Frequency [ GHz ]

BT

[ K ]

Cloudy

h=3 kmh=9 kmh=10.5 kmh=11.5 km

317 318 319 320 321 322 323 324 325

−2

0

2

4

6

8

10

Frequency [ GHz ]

∆ B

T [ K

]

Cloudy − Clear

h=3 kmh=9 kmh=10.5 kmh=11.5 km

Figure 6.7: Limb spectra at different tangent altitudes. for MASTER bandC (317 – 326 GHz). The top panel shows the clear sky spectra, the middlepanel shows the cloudy spectra and the bottom panel shows the differencebetween cloudy and clear sky spectra.

132 6 Microwave limb spectra

342 343 344 345 346 347 348 349

100

150

200

250

Frequency [ GHz ]

BT

[ K ]

Clear

h=3 kmh=9 kmh=10.5 kmh=11.5 km

342 343 344 345 346 347 348 349

100

150

200

250

Frequency [ GHz ]

BT

[ K ]

Cloudy

h=3 kmh=9 kmh=10.5 kmh=11.5 km

342 343 344 345 346 347 348 349

0

5

10

Frequency [ GHz ]

∆ B

T [ K

]

Cloudy − Clear

h=3 kmh=9 kmh=10.5 kmh=11.5 km

Figure 6.8: Limb spectra at different tangent altitudes for MASTER bandD (342 – 349 GHz). The top panel shows the clear sky spectra, the middlepanel shows the cloudy spectra and the bottom panel shows the differencebetween cloudy and clear sky spectra.

6.3 Results 133

497 498 499 500 501 502 503 504 505 506

100

150

200

250

Frequency [ GHz ]

BT

[ K ]

Clear

h=3 kmh=9 kmh=10.5 kmh=11.5 km

497 498 499 500 501 502 503 504 505 506

100

150

200

250

Frequency [ GHz ]

BT

[ K ]

Cloudy

h=3 kmh=9 kmh=10.5 kmh=11.5 km

497 498 499 500 501 502 503 504 505 506

−10

−5

0

5

10

15

Frequency [ GHz ]

∆ B

T [ K

]

Cloudy − Clear

h=3 kmh=9 kmh=10.5 kmh=11.5 km

Figure 6.9: Limb spectra at different tangent altitudes. for MASTER bandE (496 – 506 GHz). The top panel shows the clear sky spectra, the middlepanel shows the cloudy spectra and the bottom panel shows the differencebetween cloudy and clear sky spectra.

134 6 Microwave limb spectra

318 320 322 324

235

240

245

250

255

260

265

Frequency [ GHz ]

BT

[ K ]

Nadir

318 320 322 324

−1.6

−1.4

−1.2

−1

−0.8

−0.6

−0.4

−0.2

Frequency [ GHz ]

∆ B

T [ K

]

318 320 322 324

210

220

230

240

250

260

Frequency [ GHz ]

BT

[ K ]

Nadir

318 320 322 324−40

−35

−30

−25

−20

−15

−10

−5

Frequency [ GHz ]

∆ B

T [ K

]

Figure 6.10: Effect of clouds in nadir geometry. Left panel - case (1):Clear sky spectrum (grey), and cloudy spectra for Reff = 21.5 µm ( — ),Reff = 34.0 µm (– · –), Reff = 68.5 µm (· · ·), Reff = 85.5 µm (– – –) andReff = 128.5 µm ( — ). Right panel – case (4): Clear sky spectrum (grey),and cloudy spectra for IMC = 4 · 10−5 g/m3 ( — ), 1.6 · 10−3g/m3 (– · –),8 · 10−3 g/m3 (· · ·), 0.16 g/m3 (– – –) and 0.04 g/m3 ( — ) and the corre-sponding particle sizes. The top plots show absolute BTs, the bottom plotsshow the differences between cloudy and clear sky spectra.

and IMC were increased simultaneously, shows higher BT depression(see Figure 6.10).

6.4 DiscussionThe results (Figures 6.2 and 6.5) show that the ice mass content andthe size of the cloud particles both are important. Consider the dottedline in Figure 6.2, corresponding to an IMC of 1.6·10−3 g/m3 and aneffective particle radius of 64.0 µm, and the dotted line in Figure 6.5,corresponding to the same particle size but a five times higher IMC

6.4 Discussion 135

(8·10−3 g/m3): The maximum BT depression at 8 km tangent alti-tude is increased by a factor of 2.5 when we compare the two results.When the effective particle size is doubled from 34.0 to 68.5 µm, theBT depression is enhanced by a factor of 9. Consequently, in these sim-ulations the particle size is more important than the IMC. However,to estimate in detail the effect of particle size and IMC a sensitiv-ity study for more sizes and IMC is necessary. The strong impact ofboth, particle size and IMC shows, that in cloud retrievals it will bea challenge to obtain IMC and particle size simultaneously, as theyboth lead to an enhancement of the scattering effect. For retrievalswith clouds further investigations are required.

The highest effect of cloud scattering can be observed in the mi-crowave window regions of the spectral bands. The total extinctioncoefficient consists of gaseous extinction and particle extinction. Whenthe gaseous extinction cross section dominates the total extinction thescattering effect becomes smaller. In the line centers there is no differ-ence between cloudy and clear sky spectrum. At the center frequenciesthe water vapor path is so high that the transmission from the cloudto the sensor is zero, thus the existence of the cloud does not affectthe measured radiance. As gas absorption is very high at low altitudesthe cloud altitude has a big effect on the scattering signal. At low al-titudes gas absorption is dominant, therefore the scattering effect forlow clouds is very small. Cirrus clouds can exist in altitudes above10 km. Here the gas absorption is low, so the scattering signal canbe very large. The scattering coefficients increase with frequency. Butthe scattering effect does not necessarily increase with frequency. Italso depends on the gas absorption characteristics of the consideredfrequency region.

The particle shape is less important than the particle size and theIMC, at least for the cloud particles studied here. For higher aspectratios and more asymmetrical particles the impact is larger, espe-cially when the particles are oriented. For asymmetrical particles theradiation will be polarized due to scattering. First studies of the po-larization characteristics are presented in the following chapters (7and 8).

7 Simulation of polarized radiancesfor observations of cirrus cloudsin limb- and down-lookinggeometry

This chapter shows first simulations using the full capabilities of thenew scattering model. In the first part polarized simulations for a1D spherical atmosphere are presented. They show the scattering andthe polarization signal for limb- and down-looking geometries. Differ-ent particle sizes, shapes and orientations were considered. Further-more the accuracy of the scalar approximation of the radiative transferequation was validated. In the second part, polarized 3D simulationsare presented. Here it was investigated in particular, how the scat-tering and the polarization signal depend on the sensor position withrespect to the cloud. The simulations presented in this chapter havebeen published in Emde et al. (2004a).

7.1 Model simulations in a 1D sphericalatmosphere

In all simulations it was assumed that the model atmosphere consistsof nitrogen and oxygen, and the two major atmospheric trace gases:water vapor and ozone. Like in the calulations presented in the previ-ous chapter, the concentrations are taken from FASCOD (Andersonet al. (1986)) data for mid-latitudes in summer, and gas absorptionwas calculated based on the HITRAN (Rothman et al. (1998)) molecu-lar spectroscopic database using the ARTS model (version ARTS-1-0).

137

138 7 Simulation of polarized radiances

Atmospheric refraction was neglected. All calculations were carriedout for 318 GHz. The results of the calculations are summarized inTable 7.1.

7.1.1 Scattering and polarization signal fordifferent particle sizes

In order to study the impact of particle size on the radiation field,1D-calculations were carried out for four different particle sizes (equalvolume sphere radius): 25 µm, 50 µm, 75 µm and 100 µm. The parti-cles are prolate spheroids with an aspect ratio of 0.5. They are eithercompletely randomly oriented or horizontally aligned with random az-imuthal orientation. The cloud altitude is 10 – 12 km and the ice masscontent is 4.3·10−3 g/m3, which is rather small. The small value isused in order to compensate for the fact that the 1D model assumesa cloud with infinite horizontal extent. Figure 7.1 shows the radiationfield just above the cloud at 13 km altitude for completely randomlyoriented particles. The scattering signal increases significantly withthe particle size. The top panels show the difference between the scat-tered intensity field and the clearsky field. At about 90 two differentfeatures can be observed: a brightness temperature (BT) enhancementor a BT depression. The physical explanation is that the main sourceof radiation is the thermal radiation from the lower atmosphere. Forzenith angles just above 90 there is a BT enhancement because radi-ation coming from the lower atmosphere is scattered inside the cloudinto the limb directions. This part of the radiation is missing in thedown-looking directions, therefore there is a BT depression for thesedirections. The strongest scattering signal is observed in limb direc-tions, since here the path-length through the cloud is the largest. Thebottom panels of Figure 7.1 show the polarization signal, which isvery small for randomly oriented particles. The largest polarizationis observed for the largest particles (r = 100 µm) at about 91.5, buteven in this case it is below 1 K. The discrete jumps for zenith anglesfrom 100 to 180 result from the polynomial interpolation of the ra-diation field on the cloud box boundary, which is taken as radiative

7.1 Model simulations in a 1D spherical atmosphere 139

background for a clear sky calculation towards the sensor. This in-terpolation is necessary, since the intersection zenith angle of the lineof sight of the sensor with the cloud box boundary is not necessarilycontained in the optimized zenith angle grid, which is used for therepresentation of the radiation field. Since a three point polynomialinterpolation scheme is applied these jumps occur where a differentset of three points is used for the interpolation. The resolution ofthe optimized zenith angle grid is much coarser for angles close tonadir because the radiation field does not change rapidly here. Theabsolute value of the jumps is very small, they can only be seen soclearly, because the scattering signal for nadir is also very small. Theinterpolation error is below 0.2% as shown in Figure 4.7.

Figure 7.2 shows the equivalent plots for particles, which are hor-izontally aligned with random azimuthal orientation. The intensityplots are similar to the cloud case with completely randomly orientedparticles, but the polarization signal is much larger for oriented par-ticles. The maximum polarization difference (Q equals the verticalminus the horizontal intensity component) is -6.3 K for the largestparticles. In most regions Q is positive (partial vertical polarization),only in limb-directions just above 90 it is negative (partial horizontalpolarization). For randomly oriented particles, the polarization signalis due to the radiation scattered into the line of sight, because onlythe phase matrix has non-zero off-diagonal elements. For horizontallyaligned particles, the sign of the polarization signal is determined bytwo opposing mechanisms: dichroism, as manifested by a non-diagonalextinction matrix; and the effect of radiation being scattered into theline of sight. For angles just above 90 the radiation being scatteredinto the line of sight is the dominating mechanism, which results ina negative Q. For down-looking directions, where the cloud is be-tween the main radiation source and the sensor the dichroism effectis dominating, which results in a positive Q. The figure shows thatpolarization is very significant for limb radiances when the particlesare oriented.

140 7 Simulation of polarized radiances

90 95 100 105−20

−10

0

10

20

30

40

Zenith angle [ ° ]

I (cl

oudy

−cle

ar) [

K ]

25 µm50 µm75 µm100 µm

100 120 140 160 180−5

−4

−3

−2

−1

0

1

Zenith angle [ ° ]

I (cl

oudy

− c

lear

) [ K

]

25 µm50 µm75 µm100 µm

90 95 100 105−1

−0.8

−0.6

−0.4

−0.2

0

Zenith angle [ ° ]

Q [

K ]

25 µm50 µm75 µm100 µm

100 120 140 160 180−0.06

−0.05

−0.04

−0.03

−0.02

−0.01

0

0.01

Zenith angle [ ° ]

Q [

K ]

25 µm50 µm75 µm100 µm

Figure 7.1: Effect of particle size: Scattering signal of completely randomlyoriented particles with effective particle sizes 25 µm, 50 µm, 75 µm and100 µm for 318 GHz at 13 km altitude. Top panels: Intensity difference be-tween scattering calculation and clear sky calculation; bottom panels: Dif-ference between horizontal and vertical polarization.

7.1 Model simulations in a 1D spherical atmosphere 141

90 95 100 105−20

−10

0

10

20

30

40

Zenith angle [ ° ]

I (cl

oudy

−cle

ar) [

K ]

25 µm50 µm75 µm100 µm

100 120 140 160 180−5

−4

−3

−2

−1

0

1

Zenith angle [ ° ]

I (cl

oudy

− c

lear

) [ K

]

25 µm50 µm75 µm100 µm

90 95 100 105−8

−6

−4

−2

0

2

4

Zenith angle [ ° ]

Q [

K ]

25 µm50 µm75 µm100 µm

100 120 140 160 180−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Zenith angle [ ° ]

Q [

K ]

25 µm50 µm75 µm100 µm

Figure 7.2: Effect of particle size: Scattering signal of horizontally alignedparticles with effective particle sizes 25 µm, 50 µm, 75 µm and 100 µm for318 GHz at 13 km altitude. Top panels: Intensity difference between scatter-ing calculation and clear sky calculation; bottom panels: Difference betweenhorizontal and vertical polarization.

142 7 Simulation of polarized radiances

7.1.2 Effect of particle shape

In order to look at the effect of particle shape, simulations were car-ried out for particles with aspect ratios 0.5 (prolate spheroids), 1.0(spheres) and 2.0 (oblate spheroids). The particle size was 75 µmfor all calculations and ice mass content and cloud height were thesame as in the previous calculations. Figure 7.3 shows the results forcompletely randomly oriented particles. The radiation field does notchange significantly for different aspect ratios. This means that theparticle shape is not important for this particular setup. Figure 7.4shows the equivalent simulations for horizontally aligned particles withrandom azimuthal orientation. Here there are significant differencesbetween the different particle shapes. The intensity plots show thatthe BT enhancement and the BT depression are similar for all particleshapes (cf. Table 7.1). The maximum absolute values of Q are -6.4 Kand -4.0 K for oblate and prolate spheroids respectively. For sphericalparticles there is only a very small polarization signal. More simula-tions are required to study in detail the effect of particle shape on thepolarization signal.

7.1.3 Scalar simulations

In order to save CPU time and memory one can use the scalar versionof the model (cf. Section 4.1.3). To test the accuracy of the scalarapproximation, all calculations presented above were performed usingthe scalar version. Figure 7.5 shows the differences between the scalarand the vector calculations for completely randomly oriented parti-cles with different effective radii. The maximum difference for limbdirections is 0.01 K and for down-looking directions 7·10−4 K. Thesesmall differences show, that it is not necessary to use the fully polar-ized vector version to model the radiative transfer through scatteringmedia with completely randomly oriented particles. The previous sec-tion has also shown, that the polarization signal is negligible for suchcases. Figure 7.6 shows the equivalent results for horizontally alignedparticles. For down-looking directions the difference is below 0.02 K,but for limb-cases it can go up to 1.5 K. For this reason one should use

7.1 Model simulations in a 1D spherical atmosphere 143

90 95 100 105−10

−5

0

5

10

15

20

25

Zenith angle [ ° ]

I (cl

oudy

−cle

ar) [

K ]

ar 0.5ar 1.0ar 2.0

100 120 140 160 180−2.5

−2

−1.5

−1

−0.5

0

Zenith angle [ ° ]I (

clou

dy −

cle

ar) [

K ]

ar 0.5ar 1.0ar 2.0

90 95 100 105−0.7

−0.6

−0.5

−0.4

−0.3

−0.2

−0.1

0

Zenith angle [ ° ]

Q [

K ]

ar 0.5ar 1.0ar 2.0

100 120 140 160 180−0.03

−0.025

−0.02

−0.015

−0.01

−0.005

0

0.005

Zenith angle [ ° ]

Q [

K ]

ar 0.5ar 1.0ar 2.0

Figure 7.3: Effect of particle shape: Scattering signal of completely ran-domly oriented spheroidal particles with aspect ratios 0.5, 1.0 and 2.0 for318 GHz at 13 km altitude. Top panels: Intensity difference between scatter-ing calculation and clear sky calculation; bottom panels: Difference betweenhorizontal and vertical polarization.

144 7 Simulation of polarized radiances

90 95 100 105−10

−5

0

5

10

15

20

Zenith angle [ ° ]

I (cl

oudy

−cle

ar) [

K ]

ar 0.5ar 1.0ar 2.0

100 120 140 160 180−2.5

−2

−1.5

−1

−0.5

0

Zenith angle [ ° ]

I (cl

oudy

− c

lear

) [ K

]

ar 0.5ar 1.0ar 2.0

90 95 100 105−8

−6

−4

−2

0

2

4

Zenith angle [ ° ]

Q [

K ]

ar 0.5ar 1.0ar 2.0

100 120 140 160 180−0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Zenith angle [ ° ]

Q [

K ]

ar 0.5ar 1.0ar 2.0

Figure 7.4: Effect of particle shape: Scattering signal of horizontally alignedspheroidal particles with aspect ratios 0.5, 1.0 and 2.0 for 318 GHz at 13 kmaltitude. Top panels: Intensity difference between scattering calculation andclear sky calculation; bottom panels: Difference between horizontal and ver-tical polarization.

7.1 Model simulations in a 1D spherical atmosphere 145

90 95 100 105−0.01

−0.008

−0.006

−0.004

−0.002

0

Zenith angle [ ° ]

I (ve

ctor

−sca

lar)

[ K

]

25 µm50 µm75 µm100 µm

100 120 140 160 180−8

−6

−4

−2

0

2

4x 10−4

Zenith angle [ ° ]

I (ve

ctor

−sca

lar)

[ K

]

25 µm50 µm75 µm100 µm

Figure 7.5: Difference between vector RT and scalar RT calculations forcompletely randomly oriented spheroidal particles (aspect ratio 2.0) for 318GHz at 13 km altitude.

the vector model for limb RT simulations through scattering mediaconsisting of oriented particles even if one is only interested in thetotal intensity of the radiation.

90 95 100 105−1.5

−1

−0.5

0

0.5

Zenith angle [ ° ]

I (ve

ctor

−sca

lar)

[ K

]

25 µm50 µm75 µm100 µm

100 120 140 160 1800

0.002

0.004

0.006

0.008

0.01

0.012

0.014

Zenith angle [ ° ]

I (ve

ctor

−sca

lar)

[ K

]

25 µm50 µm75 µm100 µm

Figure 7.6: Difference between vector RT and scalar RT calculations forhorizontally aligned spheroidal particles (aspect ratio 2.0) for 318 GHz at13 km altitude.

146 7 Simulation of polarized radiances

Table 7.1: Summary of simulationsSize Aspect BT enh. BT dep. Polarization[µm] ratio

[−]

∆BTmax

[K]

∆BTmax

[K]

∆BT120

[K]

Qmax

[K]

Q120

[K]

p20: Completely randomly oriented particles25 7.56 -0.49 -0.04 -0.03 -0.0050 0.5 11.67 -2.61 -0.23 -0.19 -0.0075 21.32 -8.72 -0.75 -0.54 -0.01100 35.04 -19.63 -1.72 -0.96 -0.01

0.5 21.32 -8.72 -0.75 -0.54 -0.0175 1.0 20.18 -8.21 -0.70 -0.53 -0.01

2.0 21.45 -8.78 -0.76 -0.55 -0.01p30: Horizontally aligned particles with random azimuthal orientation

25 7.11 -0.47 -0.04 -1.31 0.0150 0.5 10.88 -2.47 -0.23 -2.13 0.0375 19.65 -8.18 -0.75 -4.01 0.12100 32.03 -18.29 -1.73 -6.33 0.32

0.5 19.65 -8.18 -0.75 -4.01 0.1275 1.0 19.82 -8.36 -0.71 -0.52 -0.01

2.0 18.79 -7.78 -0.75 -6.42 0.20

7.2 3D box type cloud model simulationsThe 3D version of the model was applied for simulating limb radiancesfor a cloud of finite extent embedded in a horizontally homogeneousatmosphere. The height of the cloud box was 7.3 to 12.5 km and thevertical extent of the cloud was from 9.4 to 11.5 km. The latituderange was 0 to 0.576 and the longitude range was 0 to 0.288. Thiscorresponds to a horizontal extent of approximately 64 km × 32 km.A coarse spatial discretization was chosen, because a fine resolutionis not necessary when the cloud is homogeneous; the number of gridpoints was 6 × 9 × 5. Simulations were performed for two differentIMC: 0.02 g/m3 and 0.1 g/m3 corresponding to limb optical depthsof approximately 0.5 and 2.8 respectively. The maximum propagationpath step length was set to 1 km for the optically thin cloud and to250 m for the optically thicker cloud. These values allow to assumesingle scattering for each propagation path step. It was assumed thatthe cloud consists of spheroidal ice particles with a particle size of

7.2 3D box type cloud model simulations 147

75 µm and an aspect ratio of 0.5. Calculations were performed forcompletely randomly oriented particles and for horizontally alignedparticles with azimuthally random orientation. The sensor was placedon board a satellite following a polar orbit at 820 km altitude. Ateach sensor position tangent altitudes from 0 to 13 km were measured.Figure 7.7 shows corresponding lines of sight (LOS). The figure showsthat the cloud is seen from different sides, from the top, from thebottom or from the left side. When the satellite is at a latitude of 25

the cloud is only seen for low tangent altitudes (from 0 to 6 km). Thecloud is seen at higher tangent altitudes at 27.5. For even greaterlatitudes the sensor sees the cloud from the bottom. Note that thetangent point in the first plot is behind the cloud, in the second plotin the middle of the cloud and in the last plot in front of the cloud. Inorder to compare the 1D and the 3D model versions, simulations fora 1D cloud layer with equivalent limb optical depths at 10 km tangentheight were performed. The IMC for the equivalent 1D clouds were0.005 g/m3 and 0.025 g/m3.

Figures 7.8 to 7.10 show the simulated radiances plotted as a func-tion of tangent altitude and sensor position for totally randomly ori-ented particles. The top panels show the intensity differences betweenthe clear sky calculation and the cloudy sky calculation. The bottompanels show the polarization difference Q. The contour plots on theleft hand side are the 3D results. Reddish colors indicate a brightnesstemperature enhancement due to the cloud and bluish colors indicatea BT depression. White means that there is no cloud effect. A cloudeffect can only be seen at tangent heights for which the correspondinglines of sight intersect the cloud. The intensity plot shows that up to alatitude of 27 there is a BT depression due to the cloud. The reasonis that in those cases the tangent point, from where the major sourceof thermal radiation emerges, is behind the cloud. The cloud scatterspart of the radiation away from the line of sight. For latitudes above28 a BT enhancement is observed. In these cases the tangent point isin front of the cloud. The sensor measures all radiation emerging fromthe tangent point and additionally the back-scattered radiation fromthe cloud behind the tangent point. If the tangent point is inside thecloud, between 26.5 and 28, a BT enhancement can be observed for

148 7 Simulation of polarized radiances

high tangent altitudes because part of the up-welling radiation fromthe lower atmosphere is scattered into the direction of the LOS. Forlower tangent points the scattering away from the LOS dominates,hence a BT depression is observed in this latitude range. The maxi-mum absolute values for the BT enhancement and the BT depressionare 19 K and −23 K respectively. The equivalent 1D result on the righthand side shows a larger BT enhancement of 22 K and a smaller BTdepression of −10 K. The BT depression is smaller because the op-tical depth for tangent heights below the cloud is smaller in the 1Dcalculation compared to the 3D calculation with much larger IMC.The polarization plots shows that in the 3D case as well as in the 1Dcase there is only a very small polarization difference for totally ran-domly oriented particles. In 3D, it is between −0.4 K and 0.1 K and in1D between −0.5 K and 0 K. In 1D, only negative polarization is ob-served whereas in 3D it can be positive or negative. The intensity plotin Figure 7.9 for azimuthally randomly oriented particles looks simi-lar to that for completely randomly oriented particles. However, themaximum values of BT enhancement and BT depression are slightlysmaller, about 17 K and −22 K respectively. In 1D, the intensity dif-ferences are in the range of −9 K to 20 K. The polarization differencebecomes much larger, between −3.5 K and 4.0 K can be observed inthe 3D simulation and between −4.0 K and 1.7 K in the equivalent 1Dsimulation.

Figure 7.10 shows the results of the simulation for the thicker cloudconsisting of horizontally aligned particles. The pattern looks verysimilar to that obtained for the thinner cloud but the absolute valuesof the BT depression, the BT enhancement and the polarization aremuch larger. The intensity difference is in the range from -63 K to 45 Kand the polarization difference is in the range from −7.0 K to 5.2 Kfor the 3D calculation. The equivalent 1D result ranges from −35 Kto 55 K for the intensity and from −7.0 K to 5.2 K for the polarizationdifference. Since the pattern for the thicker cloud is similar to thatobtained in the thin cloud case also for randomly oriented particles,the plot is not inlcuded here. The intensity difference ranges in thiscase from −65 K to 47 K for 3D and from −37 K to 58 K for 1D. The

7.2 3D box type cloud model simulations 149

polarization difference ranges from −0.7 K to 0.8 K for 3D and from−1.0 K to 0 K for 1D.

Overall the comparison between 1D and 3D shows similar results attangent heights inside the cloud, where the optical depth is approx-imately equivalent. For other tangent heights, the optical depths aredifferent and therefore the results deviate strongly. The scattering sig-nal in 3D depends very much on the sensor position w.r.t. the cloud.Hence it is very important to use the 3D model where the cloud extentis not very large, like in this example case, or where the clouds arehorizontally inhomogeneous.

−1010

5

10

15

20

latitude [ ° ]

heig

ht [

km ]

satellite at latitude: 25°

−1010

5

10

15

20

latitude [ ° ]

heig

ht [

km ]

satellite at latitude: 27.5°

−1010

5

10

15

20

latitude [ ° ]

heig

ht [

km ]

satellite at latitude: 29°

012345678910111213

Figure 7.7: Lines of sight (LOS) for different sensor positions and tangentaltitudes [km]. The solid line corresponds to a LOS for a tangent altitude of0 km and the dashed line to a LOS for a tangent altitude of 13 km. Dottedlines correspond to LOS for tangent heights between 0 and 13 km. Inside thesolid rectangle the single scattering properties are defined and the dashedrectangle labels the cloud box. Courtesy of Claas Teichmann.

7.2.1 Performance

The CPU time for the thin cloud cases was approximately 50 minuteson a 3 GHz Pentium 4 processor, when all four Stokes componentswere calculated. U and V are not discussed as they are approximatelyzero (less than 10−7 K) for all calculations. The computation time canbe reduced by 25% without loosing accuracy when one runs the model

150 7 Simulation of polarized radiances

Latitude [ ° ]

Tang

ent a

ltitud

e [ k

m ]

∆ I − p20

25 26 27 28 290

5

10

−20

−10

0

10

0 10 200

5

10

1D

BT [K]

Latitude [ ° ]

Tang

ent a

ltitud

e [ k

m ]

Q − p20

25 26 27 28 290

5

10

−0.3

−0.2

−0.1

0

0.1

−0.5 00

5

10

1D

BT [K]

Figure 7.8: Left panels: Scattering signal of a 3D box-type cloud embeddedin a 1D atmosphere as a function of sensor position and tangent altitudefor 318 GHz. The cloud consists of completely randomly oriented spheroidalparticles with a size of 75 µm and with an aspect ratio of 0.5. The IMC is0.02 g/m3. Right panels: 1D result for a cloud with an equivalent opticaldepth in limb (IMC = 0.005 g/m3). The upper plots show the intensity I

and the lower plots the polarization difference Q.

only for two Stokes components. The computation time for the samescenario was in this case approximately 37 minutes. The calculationfor the thicker cloud took much longer, approximately 150 minutes forall four Stokes components, because the maximum propagation-pathstep length needed to be reduced.

The computation time increases strongly with the size of the cloudbox. Doubling the number of grid points in one dimension means adoubling of the computation time. Therefore the 3D version of themodel can be used for accurate simulations to study the effect ofcloud inhomogeneity, but it is not applicable for operational use. Theperformance of the 1D version of the model is much better. All 1D

7.3 Conclusions 151

Latitude [ ° ]

Tang

ent a

ltitud

e [ k

m ]

∆ I − p30

25 26 27 28 290

5

10

−20

−10

0

10

0 10 200

5

10

1D

BT [K]

Latitude [ ° ]

Tang

ent a

ltitud

e [ k

m ]

Q − p30

25 26 27 28 290

5

10

−2

0

2

−4 −2 00

5

10

1D

BT [K]

Figure 7.9: Left panels: Scattering signal of a 3D box-type cloud embeddedin a 1D atmosphere as a function of sensor position and tangent altitudefor 318 GHz. The cloud consists of horizontally aligned spheroidal particleswith a size of 75 µm and with an aspect ratio of 0.5. The IMC is 0.02 g/m3.Right panels: 1D result for a cloud with an equivalent optical depth in limb(IMC = 0.005 g/m3). The upper plots show the intensity I and the lowerplots the polarization difference Q.

simulations shown in this chapter needed less than 30 seconds CPUtime.

7.3 ConclusionsFor the unpolarized simulations being presented in Chapter 6 as wellas for the polarized simulations presented in this chapter, the particlesize strongly influences the scattering signal. The polarization signalalso depends strongly on the particle size. Particle shape is an impor-tant cloud parameter when the cloud particles are horizontally alignedwith random azimuthal orientation. In the case of totally randomly

152 7 Simulation of polarized radiances

Latitude [ ° ]

Tang

ent a

ltitud

e [ k

m ]

∆ I − p30

25 26 27 28 290

5

10

−60

−40

−20

0

20

40

−20 0 20 400

5

10

1D

BT [K]

Latitude [ ° ]

Tang

ent a

ltitud

e [ k

m ]

Q − p30

25 26 27 28 290

5

10

−2

0

2

4

6

−5 0 50

5

10

1D

BT [K]

Figure 7.10: Left panels: Scattering signal of a 3D box-type cloud embeddedin a 1D atmosphere as a function of sensor position and tangent altitudefor 318 GHz. The cloud consists of horizontally aligned spheroidal particleswith a size of 75 µm and with an aspect ratio of 0.5. The IMC is 0.1 g/m3.Right panels: 1D result for a cloud with an equivalent optical depth in limb(IMC = 0.025 g/m3). The upper plots show the intensity I and the lowerplots the polarization difference Q.

oriented particles, changing the particle shape shows almost no effectin the simulations. For horizontally aligned particles, there is a signifi-cant difference between the scalar (unpolarized) version and the vector(polarized) version of the model in intensity. Therefore it is importantto use a vector radiative transfer model to obtain accurate results,even if one is only interested in intensity, not in polarization. The 3Dsimulations show that one must not neglect cloud inhomogeneity ef-fects. The scattering signal depends strongly upon the sensor positionwith respect to the cloud. The fact that the scattering signal is muchlarger in limb geometry compared to down-looking geometries, dueto the greater path-length through the cloud layers, demonstrates thepotential of retrieving cloud properties from limb measurements. The

7.3 Conclusions 153

major disadvantage of the 3D DOIT model is that is not yet usefulfor operational applications due to the large computation time. Butit is practical for research, for instance to study in detail the effectof different cloud parameters on polarization. The performance of the1D model is much better than that of the 3D model. Therefore the 1Dversion of the model can be applied to calculate full frequency spectraor for detailed cloud sensitivity studies.

8 A study to investigate the impactof thin layer cirrus clouds intropical regions on the EOS MLSinstrument

The Microwave Limb Sounder (MLS) experiments perform measure-ments of atmospheric composition, temperature, and pressure by limbobservations of millimeter- and sub-millimeter-wavelength thermalemission. The first MLS experiment in space was launched on theUpper Atmospheric Research Satellite (UARS) in 1991. A follow-onMLS instrument was developed for NASA’s Earth Observing System(EOS). The EOS MLS instrument was launched on the EOS Aurasatellite on the 15th of July, 2004. EOS MLS is a passive instrumentthat has radiometers in spectral bands centered near 118, 190, 240,640 and 2500 GHz. Information about the UARS and the EOS MLSinstruments is given in Waters et al. (1999). This chapter presentsscattering simulations for the EOS MLS instrument with focus onsensitivity of the scattering and polarization signal on ice mass con-tent and aspect ratio. Furthermore the possibility of retrieving shapeinformation by combination of channels with different polarizationcharacteristics is investigated.

8.1 Setup8.1.1 Selection of frequencies

The EOS MLS spectral channels, which were selected for the study,are given in Table 8.1. R1A and R1B both measure at 122 GHz, where

155

156 8 Thin layer cirrus study for EOS MLS

R1A measures the vertically polarized component of the intensity andR1B the horizontally polarized component. Hence the combination ofthe measurements of the two radiometers gives the polarization differ-ence Q at 122 GHz. Furthermore we selected the 200.5 GHz channelof R2, which measures the vertically polarized part, and the 230 GHzchannel of R3, which measures the horizontally polarized part of theradiation.

Radiometer (Polarization) Frequency Required moleculesR1A (V) 122 GHz O2, N2, H2OR1B (H) 122 GHz same as R1AR2 (V) 200.5 GHz O2, N2, H2OR3 (H) 230 GHz O2, N2, H2O, CO, O3

Table 8.1: Selected EOS MLS spectral channels for sensitivity study

8.1.2 Particle size distribution function

For all simulations the particle size distribution by Mc Farquhar andHeymsfield, which was introduced in Section 3.5.3, was used. For eachsimulation one particular particle shape was assumed. In order tostudy the effect of particle shape, many simulations were performedusing the size distribution by Mc Farquhar and Heymsfield for aspheri-cal particles with different aspect ratios. Moreover it was assumed thatthe particles are horizontally aligned.

8.1.3 Included species for absorption coefficientcalculations

Gaseous species, which have important absorption lines or continuaat the frequencies of the considered channels, are listed in Table 8.1.All of these species were included in the simulations. The selection ofgaseous species is based on Waters et al. (1999). Atmospheric profileswere taken from the FASCOD (Anderson et al., 1986) data for tropicalregions. The water vapor profile was adjusted so that the relative

8.2 Clear sky and cloudy radiances 157

humidity was 100% with respect to ice at altitudes with non-zero icemass content.

8.1.4 Definition of realistic cloud parameters

For this study box shaped clouds in pressure, latitude and longitudeare considered. The question of cloud inhomogeneity is neglected inorder to simplify the interpretation of the results with respect to theeffect of IMC, particle shape and the cloud position relative to the sen-sor. Thin layer cirrus clouds are rather homogeneous and have a largehorizontal extent, hence the homogeneous box cloud is a good approx-imation for such kind of cloud. It might even be possible to use the 1Dmodel, which is much faster than the 3D model. From observations,Del Genio et al. (2002) have derived typical cirrus cloud altitude andIMC ranges. These have been used for all simulations. The altituderange was set from 11.9 to 13.4 km. The IMC ranged from 0.0001 to1 g/m3. In order to study the effect of IMC several simulations for ho-mogeneous clouds with different IMC were performed. The horizontalextent in all 3D simulations was 400 km× 400 km, which correspondsto 3.6 latitude times 3.6 longitude in the tropics.

8.2 Clear sky and cloudy radiances at 122,200.5 and 230 GHz

As a first step the clear sky and cloudy radiances for the three channels(122, 200.5 and 230 GHz) were calculated for one special cloud caseusing the 1D model in order to see qualitatively the different behaviorof the channels. The IMC was 0.1 g/m3 and the aspect ratio was 1.5.The results are shown in Figure 8.1. The top left panel shows the clearsky radiances for tangent altitudes from 1 to 13 km. The figure shows,that for 122 GHz the gas absorption at high tangent altitudes is largercompared to 200.5 and 230 GHz so that saturation is reached at ahigher altitude. The clear sky radiances for tangent altitudes higherthan approximately 8 km are larger for 122 GHz than for 200.5 and

158 8 Thin layer cirrus study for EOS MLS

230 GHz, and for tangent altitudes below approximately 8 km the op-posite is observed. The radiances obtained for 200.5 and 230 GHz arevery similar to each other, that means that these two channels havesimilar absorption features. The top right panel shows cloudy radi-ances for the three channels and the bottom left panel shows the dif-ference between cloudy and clear sky radiances. Obviously at 122 GHzthe cloud has much less effect compared to the other channels. TheBT depression at low tangent altitudes is larger for 230 GHz com-pared to 200.5 GHz, which is mostly due to an increasing extinctioncoefficient. The bottom right panel shows the polarization differenceQ. Again this is largest for 230 GHz. Also here 200.5 and 230 GHzbehave similarly and 122 GHz looks completely different and showsa much smaller polarization signal. A positive polarization differenceis due to extinction by scattering, the radiation scattered away fromthe propagation direction is horizontally polarized, which means thatvertically polarized radiation is left in the propagation direction. Theradiation scattered into the propagation direction is horizontally po-larized. Therefore, the polarization difference Q = Iv − Ih is positive,if more radiation is scattered away from the line of sight (LOS) thaninto the LOS and negative if more radiation is scattered into the LOS.

8.3 Comparison between 1D and 3Dsimulations

In order to find out, whether it is appropriate to use the 1D model forthe thin cirrus cloud layer, calculations were performed for differentsensor positions with respect to the cloud. The LOS of the sensor atdifferent positions intersect the cloud box at different points illustratedin Figure 8.2. LOS sets A and B were also used for the comparison ofthe DOIT module with the Monte Carlo module, which was presentedin Section 5.3. The bottom panel of Figure 5.8 shows the cloud boxand the LOS sets A and B. Intersection point A is 50 km away fromthe north edge of the cloud, B is in the middle of the cloud, C is 100 kmaway from the south edge of the cloud and D is 50 km away from the

8.3 Comparison between 1D and 3D simulations 159

0 100 200 3000

5

10

15

Iclear

[K]

z tan [k

m]

100 150 200 2500

5

10

15

Icloudy

[K]

z tan [k

m]

−100 −50 0 50 1000

5

10

15

Icloudy

− Iclear

[K]

z tan [k

m]

−10 −5 0 5 100

5

10

15

Q [K]

z tan [k

m]

122 GHz200 GHz230 GHz

Figure 8.1: Top left: Clear sky radiances. Top right: Cloudy radiances.Bottom left: Cloudy minus clear sky radiances. Bottom right: Polarizationdifference.

east edge of the cloud. A single aspect ratio of 1.5 was applied for thewhole cloud. The ice mass content was assumed to be 0.1 g/m3.

The differences between the 1D and the 3D DOIT model for allLOS sets are shown in Figure 8.3. For 122 GHz the difference is forthe intensity I less than 0.2 K and for the polarization signal Q lessthan 0.1 K.

For 230 GHz and LOS set C the difference for I is more than 10 Kand for Q it is more than 4 K. For the other LOS sets the differences aremuch smaller, but still significant, especially at 12 km tangent altitude.At 13 km tangent altitude the largest difference is obtained for LOSset A and at 12 km tangent altitude the largest difference is obtained

160 8 Thin layer cirrus study for EOS MLS

Viewing direction

lat. 0.0

−1.8

−1.8lon.0.0 1.8

A (1.35,0.0)

B (0.0,0.0)

C (−0.9,0.0)

D (0.0,1.35)

1.8

Figure 8.2: Cloud box and crossing points of LOS sets A, B, C and D asseen from the top. The arrow shows the viewing direction of the instrument.

for LOS set C. This can be explained by looking at the LOS relative tothe cloud box. For A the LOS corresponding to 13 km tangent altitudeintersects the cloud latitude boundary, hence the path-length throughthe cloud is shorter in the 3D model. Therefore, at 122 GHz, wherethe cloud leads to BT depression, the 1D model yields smaller BTcompared to the 3D model. On the contrary, at 230 GHz, where thecloud leads to a BT enhancement, the 1D model yields larger BT. Thesame is valid for 12 km tangent altitude of LOS set C. Altogether,LOS set C shows larger differences than LOS set A, although theintersection point is further away from the cloud box boundary.

Apart from LOS set C the differences between the 1D and the 3Dmodel are smaller than the differences between the 3D model and the3D Monte Carlo model. Therefore, numerical inaccuracies are largerthan the error introduced by the 1D approximation. This shows, thatit is reasonable to use the 1D model for studying the effect of cloudparameters of the thin cirrus layer.

8.4 Sensitivity study 161

−0.4 −0.2 0 0.2 0.42468

101214

∆I (1D − 3D) [ K ]

z tan [

km ]

−0.05 0 0.05 0.1 0.152468

101214

∆Q (1D − 3D) [ K ]

z tan [

km ]

ABCD

Frequency: 122 GHz

Frequency: 230 GHz

Frequency: 122 GHz

Frequency: 230 GHz

−5 0 5 10 152468

101214

∆I (1D − 3D) [ K ]

z tan [

km ]

−2 0 2 4 62468

101214

∆Q (1D − 3D) [ K ]

z tan [

km ]

ABCD

Figure 8.3: Comparison of 3D simulations with 1D simulations for the dif-ferent sets of LOS. The top panels show the results for 122 GHz and thebottom panels show the results for 230 GHz.

8.4 Sensitivity study8.4.1 Dependence on ice mass content

In order to study the dependence of the total intensity and the polar-ization difference Q on the IMC, we plotted simulations for a constantaspect ratio of 3.0 and varied the IMC from 0.01 to 1 g/m3. The resultsfor all channels are presented in Figure 8.4. The left panels show thedifference between the cloudy and the clear sky intensities and theright panels show the polarization differences. For 122 GHz the BTdepression due to the cloud increases monotonically with increasingIMC. The polarization difference is positive and also increases mono-tonically with IMC. The middle and bottom panels show the resultsfor 200.5 and 230 GHz respectively. The patterns look very similarto each other but different to the 122 GHz case. For intensities, aswe have already seen in Figure 8.1, there is a BT enhancement at

162 8 Thin layer cirrus study for EOS MLS

high tangent altitudes and a BT depression at low tangent altitudes.The BT depression increases with IMC. The enhancement at 12 kmtangent altitude increases up to 0.2 g/m3, but for even higher IMCit decreases again. The reason is that the optical depth of the cloudis so large that multiple scattering events from the lower part of thecloud do not reach the top of the cloud. The polarization signal forsmall tangent altitudes is positive and increases up to approximately0.2 g/m3, then it starts decreasing. Due to multiple scattering eventsthe polarization signal is decreased. For high tangent altitudes thepolarization signal changes the sign from negative to positive at ap-proximately 0.2 g/m3. The optical depth of the cloud increases, so thatless radiation is scattered into the LOS.

8.4.2 Dependence on aspect ratio

In order to study the effect of aspect ratio on the polarized radiances,simulations for an IMC of 0.1 g/m3 are presented in Figure 8.5. Theaspect ratio is varied from 1/5 to 5. The differences between cloudyand clear sky radiances on the left hand side show that the effect ofparticle shape is very small. To study the impact of particle shapethe polarization signal is much more interesting. At all frequencies,for particles with an aspect ratio of one the polarization signal is neg-ligible and it increases in both deformation directions. At 122 GHzthe polarization signal is always positive. It increases up to approx-imately 2 K for oblate spheroids and up to approximately 1.5 K forprolate spheroids. The same behavior is seen for tangent altitudes be-low 8 km for the frequencies 200.5 and 230 GHz. However, the signal ismuch larger, up to approximately 20 K for 200.5 GHz and up to 25 Kfor 230 GHz. The absolute value of the negative polarization differenceat high tangent altitudes is smaller for 12 km tangent altitude thanbelow and above this altitude, because the path-length through thecloud is the largest in this case, which means that much radiation ismultiple scattered. This decreases the polarization signal. Again thepolarization signal is larger for oblate than for prolate particles withthe same deformation.

8.4 Sensitivity study 163

IMC [g/m3]

z tan [

km ]

I ar: 3.0 f: 122GHz

0.01 0.02 0.05 0.07 0.1 0.2 0.5 0.7 12

4

6

8

10

12

−35

−30

−25

−20

−15

IMC [g/m3]

z tan [

km ]

Q ar: 3.0 f: 122GHz

0.01 0.02 0.05 0.07 0.1 0.2 0.5 0.7 12

4

6

8

10

12

2

4

6

8

10

IMC [g/m3]

z tan [

km ]

I ar: 3.0 f: 200GHz

0.01 0.02 0.05 0.07 0.1 0.2 0.5 0.7 12

4

6

8

10

12

−100

−50

0

50

IMC [g/m3]

z tan [

km ]

Q ar: 3.0 f: 200GHz

0.01 0.02 0.05 0.07 0.1 0.2 0.5 0.7 12

4

6

8

10

12

−20

−15

−10

−5

0

5

10

IMC [g/m3]

z tan [

km ]

I ar: 3.0 f: 230GHz

0.01 0.02 0.05 0.07 0.1 0.2 0.5 0.7 12

4

6

8

10

12

−100

−50

0

50

IMC [g/m3]

z tan [

km ]

Q ar: 3.0 f: 230GHz

0.01 0.02 0.05 0.07 0.1 0.2 0.5 0.7 12

4

6

8

10

12

−20

−10

0

10

Figure 8.4: Dependence of total intensity and polarization on IMC at122 GHz, 200.5 GHz, and 230 GHz for a cloud consisting of plates with anaspect ratio of 3.0. The contours correspond to the intensity difference com-pared to clear sky radiances (left) and to the polarization difference (right).

164 8 Thin layer cirrus study for EOS MLS

aspect ratio [ ]

Tang

ent a

ltitud

e [ k

m ]

I IMC: 0.1g/m3 f: 122GHz

1/5 1/4 1/3 1/2 1 2 3 4 52

4

6

8

10

12

−3.5

−3

−2.5

−2

−1.5

−1

−0.5

aspect ratio [ ]

Tang

ent a

ltitud

e [ k

m ]

Q IMC: 0.1g/m3 f: 122GHz

1/5 1/4 1/3 1/2 1 2 3 4 52

4

6

8

10

12

0.5

1

1.5

2

aspect ratio [ ]

Tang

ent a

ltitud

e [ k

m ]

I IMC: 0.1g/m3 f: 200GHz

1/5 1/4 1/3 1/2 1 2 3 4 52

4

6

8

10

12

−20

0

20

40

60

80

aspect ratio [ ]

Tang

ent a

ltitud

e [ k

m ]

Q IMC: 0.1g/m3 f: 200GHz

1/5 1/4 1/3 1/2 1 2 3 4 52

4

6

8

10

12

−20

−10

0

10

aspect ratio [ ]

Tang

ent a

ltitud

e [ k

m ]

I IMC: 0.1g/m3 f: 230GHz

1/5 1/4 1/3 1/2 1 2 3 4 52

4

6

8

10

12

−40−20020406080

aspect ratio [ ]

Tang

ent a

ltitud

e [ k

m ]

Q IMC: 0.1g/m3 f: 230GHz

1/5 1/4 1/3 1/2 1 2 3 4 52

4

6

8

10

12

−10

0

10

20

Figure 8.5: Dependence of total intensity and polarization on aspect ratioat 122 GHz, 200.5 GHz, and 230 GHz for a cloud with a constant IMC of0.1 g/m3. The contours correspond to the intensity difference compared toclear sky radiances (left) and to the polarization difference (right).

8.5 Combination of polarized channels 165

8.5 Combination of horizontally andvertically polarized channels

Since channel R1 (122 GHz) of the EOS MLS instrument measuresboth polarizations at the same time, it might be possible to use thischannel to retrieve information about particle shape. Figure 8.6 showsscatter plots of the simulated radiances. The vertically polarized partof the radiation is plotted against the horizontally polarized part. Theresults for the particles with aspect ratio 1 (black circles) are on thediagonal, as they lead only to very small polarization. Different pointsof the same particle type correspond to different IMC. As we havealso seen in Figure 8.5 the polarization difference is positive, whichmeans that Iv is greater than Ih. The simulations show, that from themeasurements it should be possible to gain information about particleshape. The further away the measurements are from the diagonal thehigher is the deformation of the particles inside the cloud. However, itmight be difficult to distinguish between oblate and prolate particlessince both induce the same polarization state. The plot shows that forclouds with a rather large IMC the polarization difference in channel122 GHz is sufficiently large to gain information about particle shape.The scatter plot looks similar for 4 km and for 11 km tangent altitude.

Since the scattering signal is much larger in channels R2 (200.5 GHz)and R3 (230 GHz), it should also be possible to obtain informationabout particle shape from those channels, since they measure differentpolarizations. Figure 8.7 shows scattered plots of R2 (vertical polariza-tion) and R3 (horizontal polarization) at different tangent altitudes.Since the cloudy radiance is for tangent altitudes up to 9 km largerfor R2 than for R3 (compare Figure 8.1), also the points for particleswith aspect ratio 1 lie below the diagonal. The plots look similar to theplot for 122 GHz. With higher deformation the difference between R2and R3 increases. Compared to R1A/B the scattering signal is muchstronger, therefore the combination of R2 and R3 is more useful formeasuring thin cirrus clouds. For tangent altitudes inside the cloud, itdepends on the IMC, whether R2 or R3 measures larger BT. Compar-ing with Figure 8.5 it can be seen that the points above the diagonalresult from the simulations for small IMC whereas the points below

166 8 Thin layer cirrus study for EOS MLS

90 100 110 12085

90

95

100

105

110

115

120

Iv [ K ]

I h [ K

]

R1A(v) vs. R1B(h) ztan = 4 km

85 90 95 100 10585

90

95

100

105

Iv [ K ]

I h [ K

]

R1A(v) vs. R1B(h) ztan = 11 km

1/41/31/21234

Figure 8.6: Scatter plot of vertically and horizontally polarized parts ofthe intensity for 122 GHz at tangent altitudes 4 and 11 km. This resultcorresponds to measurements of channel R1A (Iv) and R1B (Ih) of theEOS MLS instrument. Different symbols correspond to cloud particles ofdifferent aspect ratios.

the diagonal result from simulations for large IMC. Again simulationsfor more extreme aspect ratios are further away from the diagonalcompared to simulations for less deformed particles. These results in-dicate that it might be possible to retrieve particle shape along withIMC from the measurements at tangent altitudes inside the cloud.

8.6 Conclusions and outlookThe first simulations for EOS MLS channels show, that the data,which will be obtained from the instrument, will be very useful tostudy cloud microphysics, like ice mass content and particle shape.Especially the different polarization characteristics can be used forthis purpose. R1A/B (122 GHz) measure both polarizations for thesame frequency at the same time, but the scattering signal is rathersmall at this frequency, so that this channel can probably only beused for studying rather thick clouds. The scattering signal in chan-nels R2 (200.5 GHz) and R3 (230 GHz) is much larger, so that thesechannels can also be used for thin clouds. Unfortunately, R2 measures

8.6 Conclusions and outlook 167

50 100 15040

60

80

100

120

140

160

Iv [ K ]

I h [ K

]R2(v) vs. R3(h) ztan = 4 km

40 60 80 100 12040

60

80

100

120

Iv [ K ]

I h [ K

]

R2(v) vs. R3(h) ztan = 8 km

50 60 70 80 9040

50

60

70

80

90

Iv [ K ]

I h [ K

]

R2(v) vs. R3(h) ztan = 9 km

50 55 60 65 7040

45

50

55

60

65

70

Iv [ K ]

I h [ K

]R2(v) vs. R3(h) ztan = 10 km

1/41/31/21234

40 50 60 7030

40

50

60

70

80

Iv [ K ]

I h [ K

]

R2(v) vs. R3(h) ztan = 11 km

20 40 60 8020

30

40

50

60

70

80

Iv [ K ]

I h [ K

]

R2(v) vs. R3(h) ztan = 12 km

Figure 8.7: Scatter plots of vertically (200.5 GHz) and horizontally(230 GHz) polarized parts of the intensity at different tangent altitudes.This result corresponds to measurements of channel R2 (200.5 GHz, Iv)and R3 (230 GHz, Ih) of the EOS MLS instrument. Different symbols cor-respond to cloud particles of different aspect ratios.

168 8 Thin layer cirrus study for EOS MLS

only vertical polarization and R3 only horizontal polarization so thatit is not possible to calculate directly the polarization difference Q foreach channel. Nevertheless it is possible to combine the two channelsin order to obtain cloud information, since absorption features andthe scattering signal in the two channels are very similar. The possi-bility in retrieving cloud information in the higher frequency channelsneeds to be investigated. The problem is these channels might be,that the atmosphere is opaque, so that the instrument can not mea-sure the clouds. The first data of the EOS MLS instrument is nowbeing analyzed in the MLS science team at the Jet Propulsion Lab-oratory (JPL). The data will be available for other science teams inthe near future. Then it will be possible to compare the DOIT modelsimulations with the real data.

9 Overall summary, conclusionsand outlook

A new scattering algorithm, called Discrete Order ITerative (DOIT)method, was developed and implemented in the Atmospheric Radia-tive Transfer Simulator (ARTS). Before starting the development, lit-erature related to radiative transfer modeling including a thermalsource and scattering was reviewed. It turned out that none of thealready existing models was well suited for the purpose of simulatingpolarized limb radiances in the microwave wavelength region, since allof the reviewed models use the plane-parallel approximation of the at-mosphere. Some models using a spherical atmosphere were developedduring the same time period as the DOIT algorithm, but these mostlyuse the scalar radiative transfer equation and other approximations.Besides the DOIT algorithm there is a Monte Carlo algorithm, whichwas also implemented into ARTS by Cory Davis in parallel to theimplementation of the DOIT algorithm. The two algorithms are atpresent the only ones, which can model polarized limb radiances.

The basic equation of the scattering model is the vector radiativetransfer equation (VRTE). The derivation of this equation requiresbasic principles and definitions of electromagnetic theory. There areseveral possibilities to calculate the scattering properties of small par-ticles. The most simple method is the Rayleigh approximation forparticles, which are very small compared to the wavelength of theradiation. For spherical particles and wavelengths comparable to theparticle size, the Lorentz-Mie theory is usually applied. Since ice par-ticles are of various shapes, which are mostly asymmetric, a moresophisticated method is required for modeling scattering of radiationin cirrus clouds. The T-matrix method, which is applicable for rota-tionally symmetric particles was chosen. Although the ice crystals are

169

170 9 Summary, conclusions and outlook

usually not rotationally symmetric, cylinders and plates are a good ap-proximation for many of the crystals. There are methods for arbitraryparticle shapes, but those methods are either computationally very ex-pensive or not well tested. The T-matrix method is currently the mostcommonly used method for the calculation of scattering properties ofcirrus clouds particles.

The VRTE is a matrix integro-differential equation which in generalcannot be solved analytically. Several numerical methods have been in-vented to solve such kind of equations; these include discrete ordinatemethods, Monte Carlo methods, and “doubling and adding” methods.The originality of the DOIT method is, that it is the first discrete or-dinate algorithm for a spherical geometry. As a platform the ARTSclear sky model was used, as it includes modules for the calculationof gas absorption, for handling ray-tracing in a three-dimensional at-mosphere and for the simulation of sensor characteristics. The DOITalgorithm solves the VRTE on a restricted part of the atmospheredenoted as the “cloud box”, in order to minimize the computationaleffort. Briefly the algorithm can be described as follows: Scattering in-tegrals, which are the difficult part of the VRTE, are first calculatedat all cloud box grid points using the clear sky field. After that theVRTE can be solved using a fixed term for the scattering integral.The solutions for all cloud box points are the first iteration radiationfield. Scattering integral fields and radiation fields are calculated alter-nately until convergence is obtained. In this way the VRTE is solvednumerically for the cloud box. The spherical geometry of the cloudbox required numerical optimizations, for instance the zenith anglegrid optimization for the representation of the radiation field.

The 1D DOIT algorithm was compared to the model FM2D devel-oped at RAL (Rutherford Appelton Laboratory) and the single scat-tering model KOPRA developed for MIPAS (Michelson Interferome-ter for Passive Atmospheric Sounding). ARTS-DOIT and the FM2Dmodel showed excellent agreement (less than 1K difference in simu-lated brightness temperatures for most cloud cases) and ARTS-DOITand KOPRA agreed well in the single scattering regime. KOPRAas well as FM2D neglect polarization. KOPRA only works for 1Dspherical atmospheres, whereas the RAL model works in 1D and 2D

171

pseudo-spherical atmospheres. The two models run faster than theARTS model, but ARTS is the more general and more accurate model.The 3D polarized DOIT algorithm was compared to the ARTS MonteCarlo algorithm. This comparison was very important, since the twoalgorithms models are the very first ones, which are able to simu-late polarized limb radiances in 3D spherical atmospheres for the mi-crowave region. The agreement between the models was satisfactory.It shows that both the Monte Carlo method and the discrete ordi-nate method can be applied for solving the VRTE in a 3D sphericalatmosphere.

Several simulation studies were performed using the new algorithm.The 1D scalar version was used to simulate frequency spectra for theMASTER instrument. The 1D polarized version was used for a sensi-tivity study of thin cirrus clouds on the EOS MLS instrument. More-over, the 3D version was used for simulations of clouds with small hor-izontal extent. The results have shown that the effect of particle sizeis very significant on both intensity and polarization of the radiation.Particle shape is an important cloud parameter when the cloud par-ticles are horizontally aligned with random azimuthal orientation. Inthe case of completely randomly oriented particles, changing the par-ticle shape shows almost no effect in the simulations. For horizontallyaligned particles, there is a significant difference between the scalar(unpolarized) version and the vector (polarized) version of the modelin intensity. Therefore it is important to use a vector radiative trans-fer model to obtain accurate results, even if one is only interested inintensity, not in polarization. The 3D simulations show that one mustnot neglect cloud inhomogeneity effects. The scattering signal dependsvery much upon the sensor position with respect to the cloud. The factthat the scattering signal is much larger in limb geometry compared todown-looking geometries, due to the greater path-length through thecloud layers, demonstrates the potential of retrieving cloud propertiesfrom limb measurements.

ARTS is a modular program and can be run in different modes. Thecomputation (CPU) time depends very much upon the chosen set-up,whether one uses the 1D- or the 3D-mode, or selects the polarized orthe unpolarized mode. CPU time can also be reduced by calculating

172 9 Summary, conclusions and outlook

two instead of all four Stokes components. The accuracy of the resultsis not affected, as long as U and V are negligible. Grid optimization isvery important for both accuracy and computation time. Overall, the1D model, with or without polarization is rather efficient and can beused for example to calculate full frequency spectra. The 3D modelhowever is very inefficient for larger 3D cloud domains, so that forsuch cases the Monte Carlo algorithm should be prefered. Althoughthe 3D calculations are computationally demanding and therefore notyet useful for operational applications, the model is practical for re-search, for instance to study in detail the effect of different cloudparameters on polarization. A feature of the DOIT method is, thatit yields the whole radiation field. To simulate radiances for differ-ent sensor positions, the radiation field only needs to be calculatedonce for the whole cloud box and the outgoing radiances can then beinterpolated on each required viewing direction.

The ARTS package, which includes besides the scattering tools(Monte Carlo and DOIT) various functions for clear sky radiativetransfer and sensor modeling, is freely available under the Gnu Gen-eral Public License and can be downloaded from http://www.sat.uni-bremen.de/arts/.

In the near future data from the EOS MLS instrument will be avail-able. This could be used first of all to validate the DOIT and MonteCarlo algorithms by testing whether they can simulate the real data.Later, after the development of a cloud retrieval algorithm, the DOITmodel can be used as a forward model for cloud parameter retrievalsfrom satellite data. The model will also be a crucial tool for the dataanalysis of the submillimeter limb sounder SMILES, which is plannedto be launched in 2008.

Appendix

A Literature review

Before starting to develop a new scattering radiative transfer modela detailed literature review about already existing models was per-formed (Emde and Sreerekha, 2004). The models that were discussedin this review are summarized in Table A.1. Atmospheric geometriesbeing used by the models and the applied methods for solving theradiative transfer equation are listed. Furthermore the table shows,which models can calculate the full Stokes vector, i.e., the polariza-tion state of scattered radiances. The last column of the table inludesthe particle types the models are able to handle.

Most of the reviewed models use a plane-parallel atmosphere. OnlySHDOM, VDISORT and the Monte Carlo model, which are three-dimensional models, use an atmosphere which is discretized using acartesian coordinate system. In these models the atmosphere does notconsist just of plane-parallel layers, but of cuboidal grid cells. Noneof the models has used a spherical geometry which is necessary tosimulate limb radiances. This is the major deficit. The new version ofARTS includes a spherical atmosphere.

Methods for 1D plane-parallel atmospheres are the Eddington ap-proximations and the doubling-and-adding method. The advantage ofthe Eddington approximations is, that they give analytical expressionsas solution and therefore the Eddington models are very fast, but theycan only be applied for spherical particles. The doubling-and-addingmethod is a simple numerical method which is also quite fast and canbe used for modeling all kinds of particle types.

The successive order of scattering method can be used for 1D and3D atmospheres, in the 3D case in combination with the discrete or-dinate method. For 3D calculations, Monte Carlo approaches are alsopossible. If one wants to calculate many viewing angles and different

175

176 A Literature review

sensor positions for a rather small scattering domain, the discrete or-dinate method is more efficient than the Monte Carlo method becausethe whole radiation field is calculated at once. On the other hand, ifonly a few viewing angles are needed, the Monte Carlo method is moreefficient.

Microwave RT models for cloudy atmospheres

Name Atmosphere Method Pol. non-sph.

MWMOD plane- DO y ySimmer (1993) parallel IterativeCzekala (1999b)RTTOV plane- D-A, n nEyre (1991) parallel EddingtonEnglish and Hewison (1998)SHDOM 3D- DO n nEvans (1998) cartesian IterativePolRadTrans plane- D-A y yEvans and Stephens (1991) parallelHybrid Model plane- Eddington n nDeeter and Evans (1998) parallel & Single

ScatteringDISORT plane- DO n nStamnes et al. (1988) parallel IterativeVDISORT plane- DO y yWeng (1992) parallel IterativeSchulz et al. (1999)Perturbation ModelGasiewski and Stalin (1990)

plane-parallel

Perturba-tionmethod

y y

Eddington ModelsKummerow (1993)

plane-parallel

Eddington n n

Appendix 177

Monte Carlo 3D- Monte y yRoberti and Kummerow(1999), Liu et al. (1996)

cartesian Carlo

VDOM 3D- DO y yHaferman et al. (1997) cartesian IterativeARTS-MC 3D- Monte y yDavis et al. (2004) spherical CarloARTS-DOIT 3D- DO y yEmde et al. (2004a) spherical Iterative

Table A.1: Overview of the reviewed radiative transfer modelsAbbreviations: Pol.– polarization, DO – Discrete-ordinate method,D-A – Doubling-and-adding method, Single Scattering – Singlescattering approximation, non-sph. – non-spherical particles

B Derivations

B.1 Solution of approximated VRTEEquation (4.7) can be solved analytically using the following matrixexponential approach

I(1) = e−〈K〉sC1 + C2, (B.1)

where C1 and C2 are constants which have to be determined. Substi-tuting (B.1) into (4.7) gives the constant C2:

−〈K〉e−〈K〉sC1 = − 〈K〉e−〈K〉sC1 − 〈K〉C2

+ 〈a〉 B +⟨S(0)

⟩C2 = 〈K〉

−1(〈a〉 B +

⟨S(0)

⟩). (B.2)

C1 can be determined using the initial condition, which is the radi-ation at the intersection point P ′ traveling towards the observationpoint P :

I(1)(s = 0) = I(0)(at intersection point) (B.3)

From the ansatz Equation (B.1) follows:

I(0) =C1 + 〈K〉−1(〈a〉 B +

⟨S(0)

⟩)C1 =I(0) − 〈K〉

−1(〈a〉 B +

⟨S(0)

⟩)(B.4)

Substituting (B.2) and (B.4) into Equation (B.1) leads to the solution:

I(1) =e−〈K〉s ·(

I(0) − 〈K〉−1(〈a〉 B +

⟨S(0)

⟩))+ 〈K〉

−1(〈a〉 B +

⟨S(0)

⟩)(B.5)

179

180 B Derivations

This can be resorted to the following form:

I(1) = e−〈K〉sI(0) +(I− e−〈K〉s

)〈K〉

−1(〈a〉 B +

⟨S(0)

⟩)(B.6)

Here I denotes the identity matrix and I(0) the Stokes vector at theintersection point. There are several ways to calculate the matrix ex-ponential functions. In ARTS the Pade-approximation is implementedaccording to Moler and Loan (1979).

B.2 Transformation of single scatteringproperties from the particle frame tothe laboratory frame for randomlyoriented particles

B.2.1 Transformation from scattering matrix tophase matrix

Instead of calculating the phase matrix Z, which relates the Stokesvectors relative to their respective meroidal planes, we can calcu-late the scattering matrix F , which relates the Stokes parametersof the incident and the scattered beams with respect to the scatter-ing plane. The scattering matrix for macroscopically isotropic andmirror-symmetric scattering media has only six independent matrixelements in contrast to the phase matrix which has in general sixteenindependent matrix elements.

From symmetry considerations follows, that the scattering matrix(Equation (3.8)) has a simple block-diagonal structure. The advan-tages of using the particle frame are obvious: On the one hand sidethe calculation of the scattering matrix using the T-matrix method isvery efficient and on the other hand side much less memory is requiredto store the phase matrix. It depends only on one angle instead of fourand it has less elements. The only draw-back is, that a transforma-tion from the particle frame to the laboratory frame is needed, as theradiative transfer calculations are performed in the laboratory frame.

B.2 Coordinate system transformation 181

Inserting the transformed ninc and nsca into Equation (3.1) we cancalculate the scattering angle

Θ = arccos(cos θsca cos θinc+sin θsca sin θinc cos(φsca−φinc)) (B.7)

where the angles θsca and φsca describe the scattered beam and θinc

and φinc the incident beam in the laboratory frame.For the transformation from the scattering matrix to the phase

matrix different cases have to be considered:1. For forward scattering (Θ = 0) the scattering frame coincides with

the laboratory frame and no transformation is required.

F = Z (B.8)

2. Different transformations are needed depending of the differencebetween azimuth angles. We can derive the following transformationformulas for 0 < φsca − φinc < π provided that θinc,sca and φinc,sca

are not equal to 0 or π:

Z(θsca, θinc, φsca, φinc) =

F11 C1F12 S1F12 0

C2F12 C1C2F22 − S1S2F33 S1C2F22 + C1S2F33 S2F34

−S1F12 −C1S2F22 − S1C2F33 −S1S2F22 + C1C2F33 C2F34

0 S2F34 −C1F34 F44

(B.9)

where

Cj = cos 2σj = 2 cos2 σj − 1 (B.10)

Sj = sin 2σj = 2√

1− cos2 σj cos σj (B.11)

j = 1, 2

The terms cos σ1 and cos σ2 can be calculated from θsca, φsca, θinc

and φinc using spherical trigonometry:

cos σ1 =cos θsca − cos θinc cos θ

sin θinc sin θ(B.12)

cos σ2 =cos θinc − cos θsca cos θ

sin θsca sin θ(B.13)

182 B Derivations

For different azimuth angles, the formulas look very similar. Onlysome signs are changed.

3. In the case that θinc,sca or φinc,sca equal zero or π, the above for-mulas are not defined. The limiting values can be derived:

limθsca→0

cos σ2 = − cos(φsca − φinc)

limθsca→π

cos σ2 = cos(φsca − φinc)

limθinc→0

cos σ1 = − cos(φsca − φinc)

limθinc→π

cos σ1 = cos(φsca − φinc) (B.14)

4. For backward scattering (Θ = π) the scattering matrix is diagonaland has only two independent elements:

F (π) =

F11(π) 0 0 0

0 F22(π) 0 00 0 −F22(π) 00 0 0 F11(π)− F22(π)

(B.15)

As the phase matrix the scattering matrix of course also depends onthe frequency and on the temperature or equivalently on the refractiveindex of the scattering medium.

B.2.2 Extinction matrix and absorption vector

For scattering media consisting of randomly oriented particles one canshow, that all off-diagonal elements of the extinction matrix K vanish.Furthermore, all diagonal elements are equal and correspond to theextinction cross-section Cext.

K = CextI = N 〈Cext〉 I (B.16)

where I is the identity matrix, N the number of particles in a volumeelement and 〈Cext〉 the average extinction cross-section per particle,which in this case is independent of the direction of propagation andof the polarization state of the incident radiation.

B.2 Coordinate system transformation 183

The ensemble-averaged emission vector for isotropic scattering me-dia must be independent of the emission direction. It can be shown,that the absorption vector just depends on the absorption cross-section Cabs

a =

Cabs

000

=

N 〈Cabs〉

000

(B.17)

where 〈Cabs〉 is the average absorption cross-section per particle.Absorption end emission cross-section depend on frequency and on

the refractive index being a function of temperature.

Acknowledgments

First of all, I would like to express my sincere thanks to my supervisorStefan Bühler for his continued interest, encouragement, support andvaluable help during the accomplishment of this work, and for goodteam work in the ARTS development. He has implemented functionsto calculate gas absorption, which are used in the DOIT algorithm.Furthermore I would like to thank my colleagues from the SAT groupfor the very nice working atmosphere and also for the useful scien-tific discussions. Special thanks goes to Sreerekha T. R. for the goodcollaboration during the development and testing phase of the DOITmodel and for the useful scientific discussions. I also appreciate thework of Claas Teichmann, who has used and tested the new modelfor his diploma thesis work on polarization. In addition I would liketo thank Oliver Lemke for programming assistance and technical sup-port.

Next, I would like to thank Cory Davis from the University of Ed-inburgh, who has developed the Monte Carlo scattering model. Sincemany parts of the models are in common, our collaboration was veryexpedient. In addition, I would like to thank Patrick Eriksson fromthe Chalmers Institute of Technology in Gothenburg, who is also partof the ARTS development team, for suggestions and ideas. He hasimplemented the 3D ray-tracing scheme, which is used in the DOITalgorithm.

A large fraction of the work was carried out in the context of anESTEC study (Contract No. 15457/02/NL/MW). During this study,the comparisons of the DOIT model with the models FM2D and KO-PRA were carried out. I appreciate very much the collaboration withRichard Siddans from the the Rutherford Appleton Laboratory (RAL)near Oxford, who has implemented scattering in the model FM2D

185

186

and who helped me finding several bugs in the initial DOIT algo-rithm. Furthermore, I would like to thank Michael Höpfner from the“Forschungszentrum” in Karlsruhe, who did a lot of work for the com-parison study between the DOIT model and the KOPRA model, inwhich he has included a scattering algorithm. I also thank all othermembers of the consortium for the good team work. I very much ap-preciate suggestions and comments from Brian Kerridge (RAL), whowas the supervisor of the science consortium, and from Jörg Langenfrom ESTEC, whose criticism was always well founded and helpful.

For the calculation of single scattering properties, many differentpublic domain programs were used. I would like to thank MichaelMishchenko and Steven Warren for making available the T-matrixprogram and the refractive index program, respectively. FurthermoreI thank Christian Mätzler for providing the Mie program.

I would like to express my gratitude to Prof. Klaus Künzi and Prof.Clemens Simmer for reviewing this thesis and for helpful commentsand suggestions.

Last but not least I would like to thank my parents and my friends,who were always there when I needed them. Special thanks goes toBruno Matzas for reading and commenting on the manuscript.

Besides the ESTEC study mentioned earlier, this work was fundedby the German Federal Ministry of Education and Research (BMBF),within the DLR project SMILES, grant 50 EE 9815, and within theAFO2000 project UTH-MOS, grant 07ATC04. It is also a contributionto the COST Action 723 ‘Data Exploitation and Modeling for theUpper Troposphere and Lower Stratosphere’.

C List of acronyms

Acronym MeaningARTS Atmospheric Radiative Transfer SimulatorBT Brightness TemperatureCLEAS Cryogenic Limb Array Etalon SpectrometerCPU Central Processing UnitCRISTA Cryogenic Spectrometers and Telescopes for the

AtmosphereDDA Discrete Dipole ApproximationDDSCAT Discrete Dipole Approximation for Scattering

and Absorption of Light by Irregular ParticlesDISORT Discrete Ordinate Radiative Transfer ModelDOIT Discrete Ordinate ITerative methodDOM Discrete Ordinate MethodECBM Extended Boundary Condition MethodEOS Earth Observing SystemESA European Space AgencyESTEC European Space research and TEchnology Cen-

treFASCOD Fast Atmosphere Signature CodeFIRE First ISCCP Regional ExperimentFM2D Forward Model 2DGCM Global Climate ModelsGOME Global Ozone Monitoring ExperimentGOMETRAN GOME radiative TRANsfer modelHITRAN High-resolution Transmission Molecular Absorp-

tion database

187

188 C List of acronyms

Acronym MeaningISCCP International Satellite Cloud Climatology

ProjectIMC Ice Mass ContentJPL Jet Propulsion LaboratoryKOPRA Karlsruhe Optimized and Precise Radiative

Transfer AlgorithmLOS Line Of SightMATLAB MAtrix LABoratoryMAS Millimeter Atmospheric SounderMASTER Millimeter Wave Acquisitions for Strato-

sphere/Troposphere Exchange ResearchMC Monte Carlo methodMIPAS Michelson Interferometer for Passive Atmo-

spheric SoundingMLS Microwave Limb SounderMWMOD MicroWave MODelNASA National Aeronautics and Space AdmisistrationPyARTS Python ARTSRAL Rutherford Appleton LaboratoryRT Radiative TransferRTTOV fast Radiative Transfer model for TOVsSHDOM Spherical Harmonics Discrete Ordinate MethodSMILES Superconduction Submillimeter-Wave Limb

Emission SounderSMR Submillimeter RadiometerSRTE Scalar Radiative Transfer EquationSWCIR Submillimeter-Wave Cloud Ice RadiometerTES Troposhperic Emission SpectrometerTOVS Tiros Operational Vertical SounderUARS Upper Atmospheric Research SatelliteVDISORT Vector Discrete Ordinate Radiative Transfer

ModelVDOM Vector Discrete-Ordinates MethodVRTE Vector Radiative Transfer Equation

D List of symbols

Symbol Definition and dimension in SI units Introducedinsection

ap particle absorption vector [m2] 1.3.4〈a〉 total ensemble averaged absorption vector

(includes particle and gas contributions)[m−1]

1.5

〈ag〉 averaged gas absorption vector [m−1] 1.5〈ap〉 ensemble averaged particle absorption

vector [m−1]1.5

〈api 〉 ensemble averaged absorption vector for

one particle type [m2]1.5

B Planck blackbody energy distribution[W s m−2 sr−1]

1.3.4

c speed of light in vacuum [m s−1] 1.1Cabs absorption cross section [m2] 1.3.5Cext extinction cross section [m2] 1.3.5Csca scattering cross section [m2] 1.3.5E electric field vector [V m−1] 1.1E0 amplitude of electric field vector [V m−1] 1.1Eθ, Eφ spherical coordinate components of the

electric field vector [V m−1]1.2

Eh, Ev horizontal and vertical components of theelectric field vector [V m−1]

1.2

F scattering matrix [m2] 3.3g probability density function [–] 5.3H magnetic field vector [A m−1] 1.1

189

190 D List of symbols

Symbol Definition and dimension in SI units Introducedinsection

~ Planck constant divided by 2π [Js] 1.3.4I intensity, first Stokes parameter [W m−1] 1.2I specific intensity [W s m−2 sr−1] 1.5I Stokes vector [W m−1] 1.2I specific intensity vector or “Stokes vector”

[W s m−2 sr−1]1.5

Ib blackbody Stokes column vector[W s m−2 sr−1]

1.3.4

I radiation field [W s m−2 sr−1] 4.1.1I(n) nth order radiation field [W s m−2 sr−1] 4.1.1IMC ice mass content [kg m−3] 3.5k = kR +kI

(complex) wave number [m−1] 1.1

kb Boltzmann constant [JK−1 ] 1.3.4K total extinction matrix [m−1] 1.3.3〈K〉 ensemble averaged total extinction matrix

[m−1]1.5

〈Kg〉 ensemble averaged gaseous extinction ma-trix [m−1]

1.5

〈Kp〉 ensemble averaged particle extinction ma-trix [m−1]

1.5

〈Kpi 〉 ensemble averaged extinction matrix for

one particle type [m2]1.5

m =mR + mI

(complex) refractive index relative to vac-uum of surrounding medium [–]

1.1

m mass of a particle [kg] 3.5N number of particles [–] 1.4ng volume mixing ratio [–] 1.5np particle number density [m−3] 1.5n(r) particle size distribution function [–] 3.5n unit vector [–] 1.2ninc unit vector in the incidence direction [–] 1.3

191

Symbol Definition and dimension in SI units Introducedinsection

nsca unit vector in the scattering direction [–] 1.3p phase function [–] 1.3.5p degree of polarization [–] 1.2plin degree of linear polarization [–] 1.2pcirc degree of circular polarization [–] 1.2P time-averaged Poynting vector [W m−2] 1.1~p pressure grid [Pa] 4.1.1Q second Stokes parameter [W m−1] 1.2Q second component of specific intensity

vector[W s m−2 sr−1]

1.5

Qabs absorption efficiency [–] 3.3Qext extinction efficiency [–] 3.3Qsca scattering efficiency [–] 3.3r distance from the origin of a coordinate

system [m]1.3

r equal volume sphere radius of a particle[m]

3.2

r random number [–] 5.3r radius (position) vector [m] 1.1Reff effective radius of a particle size distribu-

tion [m]3.5

Rme median radius of a particle size distribu-tion

5.1

〈S〉 scattering source function [W s m−1 sr−1] 5.1S(n) nth order scattering integral field [W s

m−2 sr−1]4.1.1

t time [s] 1.1T temperature [K] 1.3.4TPlanck Planck brightness temperature [K] 2.7TRJ Rayleigh Jeans brightness temperature

[K]2.7

192 D List of symbols

Symbol Definition and dimension in SI units Introducedinsection

U third Stokes parameter [W m−1] 1.2U third component of specific intensity vec-

tor[W s m−2 sr−1]

1.5

V fourth Stokes parameter [W m−1] 1.2V fourth component of specific intensity vec-

tor[W s m−2 sr−1]

1.5

V volume of a particle [m−3] 3.5W power [W] 1.3.4x scattering parameter [–] 3.3Z phase matrix [m2] 1.3.2〈Z〉 ensemble averaged phase matrix [m−1] 1.4〈Zi〉 ensemble averaged phase matrix for one

particle type [m2]1.5

α polarizability [m3] 3.3αg

i individual gas absorption coefficient [m−1] 1.5αp particle absorption coefficient [m−1] 1.1〈αg〉 averaged gas absorption coefficient [m−1] 1.5~α latitude grid [] 4.1.1~β longitude grid [] 4.1.1∆s path length element [m] 5.1∆S surface element [m2] 1.3.3∆ω angular frequency interval [s−1] 1.3.4∆Ω solid angle [sr] 1.3.4∆n phase of amplitude matrix [–] 1.4ε electric permittivity [F m−1] 1.1λ free space wavelength [m] 1.1µ magnetic permeability [H m−1] 1.1ν frequency of radiation [s−1] 1.1ω angular frequency [s−1] 1.1ω0 single scattering albedo [–] 1.3.5

193

Symbol Definition and dimension in SI units Introducedinsection

ω similar to scattering albedo, used inMonte Carlo model [–]

5.3

π pi [–] 1.5τmax maximal optical depth [–] 4.1.4~φ azimuth angle grid [] 4.1.1ρ density of a scattering medium [gm−3] 3.5~θ zenith angle grid [] 4.1.1Θ scattering angle [] 3.2

194 D List of symbols

Symbol Definition and dimension in SI units Introducedinsection

General notation

x∗ complex conjugate of x 1.1〈x〉 ensemble average of x 1.4X matrix X 1.1Xij element (ij) of X 1.3.2x vector x 1.1xi ith element of x 1.5xinc x for incident direction 1.3.1xsca x for scattering direction 1.3.1XT transpose of X 1.5∫4π

integral over the whole space 1.5Rex Real part of x 1.1Imx Imaginary part of x 1.1x spatial average of x 2.3.3Γ gamma function 3.5I identity matrix B.1

Bibliography

Anderson, G. P., Clough, S. A., Kneizys, F. X., Chetwynd, J. H.and Shettle, E. P., 1986: AFGL atmospheric constituent profiles(0–120 km). Tech. Rep. TR-86-0110, AFGL.

Arking, A., 1991: The radiative effects of clouds and their impact onclimate. Bull. Am. Meteorol. Soc., 72, 795–813.

Beer, R., Glavich, T. A. and Rider, D. M., 2001: Tropospheric emissionspectrometer for the Earth Observing System’s Aura satellite. Appl.Opt., 40, 15, 2356–2367.

Bohren, C. and Huffman, D. R., 1998: Absorption and Scattering ofLight by Small Particles. Wiley Science Paperback Series.

Buehler, S. A., 1999: Microwave Limb Sounding of the Stratosphereand Upper Troposphere. Berichte aus der Physik, PO-Box 1290, D52013 Aachen: Shaker Verlag GmbH. PhD thesis, University of Bre-men, ISBN 3-8265-4745-4, 262 pages.

Buehler, S. A., Eriksson, P., Kuhn, T., von Engeln, A. and Verdes,C., 2005a: ARTS, the Atmospheric Radiative Transfer Simulator.J. Quant. Spectrosc. Radiat. Transfer, 91, 1, 65–93.

Buehler, S. A., Verdes, C. L., Tsujimaru, S., Kleinboehl, A., Bremer,H., Sinnhuber, M. and Eriksson, P., 2005b: The Expected Perfor-mance of the SMILES Submillimeter-Wave Limb Sounder comparedto Aircraft Data. Radio Sci. In press.

von Clarmann, T., Dudhia, A., Edwards, D. P., Funke, B., Höpfner,M., Kerridge, B., Kostsov, V., Linden, A., López-Puertas, M. andTimofeyev, Y. M., 2002: Intercomparison of radiative transfer codesunder non-local thermodynamic equilibrium conditions. J. Geophys.Res., 107, D22, 4631, doi:10.1029/2001JD001 551.

von Clarmann, T., Höpfner, M., Funke, B., López-Puertas, M., Dud-hia, A., Jay, V., Schreier, F., Ridolfi, M., Ceccherini, S., Kerridge,

195

196 BIBLIOGRAPHY

B. J., Reburn, J. and Siddans, R., 2003: Modelling of atmosphericmid–infrared radiative transfer: the AMIL2DA algorithm intercom-parison experiment. J. Quant. Spectrosc. Radiat. Transfer, 78, 3-4,381–407 doi:10.1016/S0022–4073(02)00 262–5.

Czekala, H., 1999a: Microwave radiative transfer calculations withmultiple scattering by nonspherical hydrometeors. Ph.D. thesis,Rheinische Friedrich-Wilhems-Universität Bonn, Auf dem Hügel 20,53121 Bonn.

Czekala, H., 1999b: Microwave radiative transfer calculations withmultiple scattering by nonspherical hydrometeors. Ph.D. thesis,Rheinische Friedrich-Wilhelms-Universitaet Bonn, Auf dem Huegel20, 53121 Bonn, Germany.

Czekala, H. and Simmer, C., 1998: Microwave radiative transfer withnonspherical precipitating hydrometeors. J. Quant. Spectrosc. Ra-diat. Transfer, 60, 3, 365–374.

Davis, C., Emde, C. and Harwood, R., 2004: A 3D Polarized ReversedMonte Carlo Radiative Transfer Model for mm and sub-mm PassiveRemote Sensing in Cloudy Atmospheres. IEEE T. Geosci. Remote.In press.

Deeter, M. N. and Evans, K. F., 1998: A hybrid Eddington-single scat-tering radiative transfer model for computing radiances from ther-mally emitting atmospheres. J. Quant. Spectrosc. Radiat. Transfer,60, 4, 635–648.

Del Genio, A. D., Wolf, A. B., Mace, G. G. and Miloshevich, L. M.,2002: Observed Regimes of Mid-latitude and Tropical Cirrus Mir-cophysical Behavior. In Twelfth ARM Science Team Meeting Pro-ceedings, St. Petersburg, Florida.

Donovan, D. P., 2003: Ice-cloud effective particle size parametriza-tion based on combined lidar, radar reflectivity, and meanDoppler velocity measurements. J. Geophys. Res., 108, D18, 4573.Doi:10.1029/2002JD003469.

Draine, B. T. and Flateau, P. J., 2003: User Guide to the DiscreteDipole Approximation Code DDSCAT.6.0. url http://arxiv.org/abs/astro-ph/0300969.

Emde, C., Buehler, S. A., Davis, C., Eriksson, P., Sreerekha, T. R. andTeichmann, C., 2004a: A Polarized Discrete Ordinate Scattering

BIBLIOGRAPHY 197

Model for Simulations of Limb and Nadir Longwave Measurementsin 1D/3D Spherical Atmospheres. J. Geophys. Res., 109, D24.

Emde, C., Buehler, S. A., Eriksson, P. and Sreerekha, T. R., 2004b:The effect of cirrus clouds on microwave limb radiances. J. Atmos.Res., 72, 1–4, 383–401.

Emde, C. and Sreerekha, T. R., 2004: Development of a RT modelfor frequencies between 200 and 1000 GHz, WP1.2 Model Review.Tech. Rep., ESTEC Contract No AO/1-4320/03/NL/FF.

English, S. and Hewison, T., 1998: A fast generic millimeter waveemissivity model. In Microwave Remote Sensing of the Atmosphereand the Environment, T. Hayasaka, Y. J., D.L. Wu and Jiang, J.,eds., vol. 3503 of Proceedings of SPIE, SPIE.

Eriksson, P. and Buehler, S. A. (Eds.), 2001: Atmospheric Millimeterand Sub-Millimeter Wave Radiative Transfer Modeling II. LogosVerlag Berlin. ISBN 3-89722-585-9, ISBN 1615-6862.

Eriksson, P., Buehler, S. A., Emde, C., Sreerekha, T. R., Melsheimer,C. and Lemke, O., 2004: ARTS-1-1 User Guide. Universityof Bremen. 308 pages, regularly updated versions available atwww.sat.uni-bremen.de/arts/.

Eriksson, P., Merino, F., Murtagh, D., Baron, P., Ricaud, P. and de laNoë, J., 2002: Studies for the Odin sub-millimetre radiometer: 1.Radiative transfer and instrument simulation. Can. J. Phys., 80,321–340.

Evans, K. F., 1998: The spherical harmonics disrete ordinate methodfor three-dimensional atmospheric radiative transfer. J. Atmos. Sci.,55, 429–466.

Evans, K. F. and Stephens, G. L., 1991: A new polarized atmosphericradiative transfer model. J. Quant. Spectrosc. Radiat. Transfer, 46,5, 412–423.

Evans, K. F., Walter, S. J., Heymsfield, A. J. and Deeter, M. N.,1998: Modeling of Submillimeter Passive Remote Sensing of CirrusClouds. J. Appl. Met., 37, 184–205.

Evans, K. F., Walter, S. J., Heymsfield, A. J. and McFarquhar,G. M., 2002: Submillimeter-Wave Cloud Ice Radiometer: Simula-tions of retrieval algorithm performance. J. Geophys. Res., 107,D3, 10.1029/2001JD000 709.

198 BIBLIOGRAPHY

Eyre, J. R., 1991: A fast radiative transfer model for satellite soundingsystems. Tech. Rep., ECMWF Tech. Memo. No. 176.

Fischer, H. and Oelhaf, H., 1996: Remote sensing of vertical profilesof atmospheric trace constituents with MIPAS limb-emission spec-trometers. Appl. Opt., 35, 16, 2787–2796.

Gasiewski, A. J. and Stalin, D. H., 1990: Numerical modeling of pas-sive microwave O2 observations over precipitation. Radio Science,25, 3, 217–235.

Glatthor, N., Höpfner, M., Stiller, G. P., von Clarmann, T., Dud-hia, A., Echle, G., Funke, B. and Hase, F., 1999: Intercomparisonof the KOPRA and the RFM Radiative Transfer Codes. In Proc.European Symposium on Atmospheric Measurements from Space,ESAMS’99, 18–22 Jan 1999, Noordwijk, European Space Agency,ESTEC, Noordwijk, The Netherlands.

Haferman, J. L., Smith, T. F. and Krajewski, W. F., 1997: A multi-dimensional discrete-ordinates method for polarised radiative trans-fer. Part 1: Validation for randomly oriented axisymmetric particles.J. Quant. Spectrosc. Radiat. Transfer, 58, 379–398.

Hartmann, G. K., Bevilacqua, R. M., Schwartz, P. R., Kämpfer, N.,Kuenzi, K. F., Aellig, C. P., Berg, A., Boogaerts, W., Connor,B. J., Croskey, C. L., Daehler, M., Degenhardt, W., Dicken, H. D.,Goldizen, D., Kriebel, D., Langen, J., Loidl, A., , Olivero, J. J.,Pauls, T. A., Puliafito, S. E., Richards, M. L., Rudin, C., Tsou,J. J., Waltman, W. B., Umlauft, G. and Zwick, R., 1996: Measure-ments of O3 , H2O, and ClO in the middle atmosphere using theMillimeter–wave Atmospheric Sounder (MAS). Geophys. Res. Lett.,23, 17, 2313–2316.

Havemann, S. and Baran, A. J., 2001: Extension of T-matrix to scat-tering of electromagnetic plane waves by non-axisymmetric dielec-tric particles: application to hexagonal ice cylinders. J. Quant. Spec-trosc. Radiat. Transfer, 70, 2, 139–158.

Heymsfield, A. J. and Platt, C. M. R., 1984: A parametrization ofthe particle size spectrum of ice clouds in terms of the ambienttemperature and the ice water content. J. Atmos. Sci., 41, 864–855.

Hoepfner, M. and Emde, C., 2005: Comparison of single and multiple

BIBLIOGRAPHY 199

scattering approaches for the simulation of limb-emission observa-tions in the mid-IR. J. Quant. Spectrosc. Radiat. Transfer, 91, 3,275–285.

Ishimoto, H. and Masuda, K., 2002: A Monte Carlo approach forthe calculation of polarized light: application to an incident nar-row beam. J. Quant. Spectrosc. Radiat. Transfer, 72, 4, 462–483.

Jackson, J. D., 1998: Classical electrodynamics. New York: John Wiley& Sons.

Kerridge, B., Jay, V., Reburn, J., Siddans, R., Latter, B., Lama,F., Dudhia, A., Grainger, D., Burgess, A., Höpfner, M., Steck, T.,Stiller, G., Buehler, S., Emde, C., Eriksson, P., Ekström, M., Baran,A. and Wickett, M., 2004: Consideration of mission studying chem-istry of the UTLS, Task 2 Report. Tech. Rep., ESTEC Contract No15457/01/NL/MM.

Kinne, S. et al., 1997: Cirrus cloud radiative and microphysical prop-erties from ground observations and in situ measurements duringFIRE 1991 and their applications to exhibit problems in cirrus so-lar radiative transfer modeling. Journal of Atmospheric Science, 54,2320–2344.

Kummerow, C., 1993: On the accuracy of the Eddington Approxima-tion for radiative transfer in the microwave frequencies. J. Geophys.Res., 98, D2, 2757–2765.

Liou, K. N., 2002: An Introduction to Atmospheric Radiation, vol. 84of International Geophysics Series. Academic Press, Elsevier Sci-ence.

Liou, K. N., Yang, P., Takano, Y., Sassen, K., Charlock, T. P. andArnott, W. P., 1998: On the radiative properties of contrail cirrus.Geophys. Res. Lett., 25, 1161-1164.

Liu, Q., Simmer, C. and Ruprecht, E., 1996: Three-dimensional ra-diative transfer effects of clouds in the microwave spectral range. J.Geophys. Res., 101, D2, 4289–4298.

Mätzler, C., 1998: Solar System Ices, vol. 227 of Astrophys. and SpaceSci. Library, chap. Microwave properties of ice and snow. Dordrecht:Kluwer Academic publishers, 241–257.

Mätzler, C., 2002: MATLAB Functions for Mie Scattering and Ab-sorption - Version 2. Tech. Rep. 2002-11, Universität Bern.

200 BIBLIOGRAPHY

McFarquhar, G. M. and Heymsfield, A. J., 1997: Parametritation ofTropical Cirrus Ice Crystal Size Distributions and Implications forRadiative Transfer: Results from CEPEX. J. Atmos. Sci., 54, 2187–2200.

Mishchenko, M. I., 2000: Calculation of the amplitude matrix fora nonspherical particle in a fixed orientation. Applied Optics, 39,1026–1031.

Mishchenko, M. I., 2002: Vector radiative transfer equation for ar-bitrarily shaped and arbitrarily oriented particles: a microphysicalderivation from statistical electromagnetics. Appl. Opt., 41, 7114–7134.

Mishchenko, M. I., Hovenier, J. W. and Travis, L. D. (Eds.), 2000:Light Scattering by Nonspherical Particles. Academic Press, ISBN0-12-498660-9.

Mishchenko, M. I., Travis, L. D. and Lacis, A. A., 2002: Scattering,Absorption, and Emission of Light by Small Particles. CambridgeUniversity Press.

Moler, C. B. and Loan, C. F. V., 1979: Nineteen Dubious Ways toCompute the Exponential of a Matrix. SIAM Review, 20, 801–836.

Murtagh, D., Frisk, U., Merino, F., Ridal, M., Jonsson, A., Stegman,J., Witt, G., Eriksson, P., Jimenez, C., Megie, G., de la Noe, J.,Ricaud, P., Baron, P., Pardo, J. R., Hauchcorne, A., Llewellyn,E. J., Degenstein, D. A., Gattinger, R. L., Lloyd, N. D., Evans,W. F. J., McDade, I. C., Haley, C. S., Sioris, C., von Savigny,C., Solheim, B. H., McConnell, J. C., Strong, K., Richardson,E. H., Leppelmeier, G. W., Kyro, E., Auvinen, B. H. and Oikari-nen, L., 2002: An overview of the Odin atmospheric mission. Can.J. Phys., 80, 309–319, url www.nrc.ca/cgi-bin/cisti/journals/rp/rp2_abst_e?cjp_p01-157_80_ns_nf_cjpS1-02.

Oikarinen, L., Sihvola, E. and Kyrola, E., 1999: Multiple scatteringradiance in limb-viewing geometry. J. Geophys. Res., 104, D24,31 261–31 274.

Roberti, L., Haferman, J. and Kummerow, C., 1994: Microwave ra-diative transfer through horizontally inhomogeneous precipitatingclouds. J. Geophys. Res., 99, D8, 16,707–16,718.

Roberti, L. and Kummerow, C., 1999: Monte Carlo calculations of

BIBLIOGRAPHY 201

polarized microwave radiation emerging from cloud structures. J.Geophys. Res., 104, D2, 2093–2104.

Roche, A. E., Kumer, J. B., Mergenthaler, J. L., Ely, G. A., Uplinger,W. G., Potter, J. F., James, T. C. and Sterritt, L. W., 1993: TheCryogenic Limb Array Etalon Spectrometer CLAES on UARS: Ex-periment Description and Performance. J. Geophys. Res., 98, D6,10,763–10,775.

Rothman, L. S., Rinsland, C. P., Goldman, A., Massie, S. T., Ed-wards, D. P., Flaud, J.-M., Perrin, A., Camy-Peyret, C., Dana,V., Mandin, J.-Y., Schroeder, J., McCann, A., Gamache, R. R.,Wattson, R. B., Yoshino, K., Chance, K. V., Jucks, K. W., Brown,L. R., Nemtchinov, V. and Varanasi, P., 1998: The HITRAN Molec-ular Spectroscopic Database and HAWKS (HITRAN AtmosphericWorkstation): 1996 edition. J. Quant. Spectrosc. Radiat. Transfer,60, 665–710.

Rozanov, V. V., Diebel, D., Spurr, R. J. D. and Burrows, J. P., 1997:GOMETRAN: A radiative transfer model for the satellite projectGOME - the plane-parallel version. . J. Geophys. Res. 102 (d14) ,16683-16695.

Schulz, F. M. and Stamnes, K., 2000: Angular distribution of theStokes vector in a plane-parallel, verically inhomogeneous mediumin the vector discrete ordinate radiative transfer. J. Quant. Spec-trosc. Radiat. Transfer, 65, 609–620.

Schulz, F. M., Stamnes, K. and Wenig, F., 1999: An improved and gen-eralized Discrete Ordinate Method for polarized (vector) radiativetransfer. J. Quant. Spectrosc. Radiat. Transfer, 61, 1, 105–122.

Simmer, C., 1993: Handbuch zu MWMOD, Kieler Strahlungstrans-portmodell für den Mikrowellenbereich. Institut für MeereskundeKiel.

Skofronik-Jackson, G. M., Gasiewski, A. J. and Wang, J. R., 2002:Influence of Microphysical Cloud Parametrizations on Microwavebrightness Temperatures. IEEE Transactions on Geoscience andRemote Sensing, 40, 1.

Spang, R., Riese, M., Eidmann, G., Offermann, D. and Wang, P. H.,2001: A detection method for cirrus clouds using CRISTA 1 and 2measurements. Adv. Space Res., 27, 10, 1629–1634.

202 BIBLIOGRAPHY

Sreerekha, T. R., Buehler, S. A. and Emde, C., 2002: A Simple NewRadiative Transfer Model for Simulating the Effect of Cirrus Cloudsin the Microwave Spectral Region. J. Quant. Spectrosc. Radiat.Transfer, 75, 611–624.

Stallman, R. M., 2002: Free Software, Free Society: Selected Essays.GNU Press, ISBN 1-882114-98-1.

Stamnes, K., Tsay, S.-C., Wiscombe, W. and Jayaweera, K., 1988:Numerically stable algorithm for discrete-ordinate-method radiativetransfer in multiple scattering and emitting layered media. Appl.Opt., 27, 2502–2509.

Stiller, G. P. (Ed.), 2000: The Karlsruhe Optimized and Precise Ra-diative Transfer Algorithm (KOPRA), vol. FZKA 6487 of Wis-senschaftliche Berichte. Forschungszentrum Karlsruhe.

Stiller, G. P., von Clarmann, T., Funke, B., Glatthor, N., Hase, F.,Höpfner, M. and Linden, A., 2002: Sensitivity of trace gas abun-dances retrievals from infrared limb emission spectra to simplifyingapproximations in radiative transfer modeling. J. Quant. Spectrosc.Radiat. Transfer, 72, 249–280.

Stubenrauch, C. J., Holz, R., Chedin, A., Mitchell, D. L. and Baran,A. J., 1999: Retrieval of cirrus ice crystal sizes from 8.3 and 11.1 mmemissivities determined by the improbed initialization inversion ofTIROS-N Operational Vertical Sounder observations. J. Geophys.Res., 104, D24, 31 793–31 808.

Tjemkes, S. A., Patterson, T., Rizzi, R., Shephard, M. W., Clough,S. A., Matricardi, M., Haigh, J. D., Höpfner, M., Payan, S., Trot-senko, A., Scott, N., Rayer, P., Taylor, J. P., Clerbaux, C., Strow,L. L., DeSouza-Machado, S., Tobin, D. and Knuteson, R., 2003:The ISSWG line–by–line inter–comparison experiment. J. Quant.Spectrosc. Radiat. Transfer, 77, 4, 433–453.

Toon, O. B., Tolbert, M. A., Middlebrook, A. M. and Jordan, J., 1994:Infrared optical constants of H2O, ice, amorphous nitric acid solu-tions, and nitric acid hydrates. J. Geophys. Res., 99, D12, 25 631–25 654.

van de Hulst, H. C., 1957: Light Scattering by Small particles. NewYork: Dover Publications. Corrected republication 1981.

Wang, P.-H., Minnis, P., McCormick, M. P., Kent, G. S. and Skeens,

BIBLIOGRAPHY 203

K. M., 1996: A 6-year climatology of cloud occurence frequency fromSratospheric Aerosol and Gas Experiment II observations (1985–1990). J. Geophys. Res., 101, 29, 407–429.

Warren, S., 1984: Optical Constants of Ice from the Ultraviolet to theMicrowave. Appl. Opt., 23, 1206–1225.

Waterman, P. C., 1965: Matrix formulation of electromagnetic scat-tering. Proc. IEEE, 53, 805–812.

Waters, J. W., Read, W. G., Froidevaux, L., Jarnot, R. F., Cofield,R. E., Flower, D. A., Lau, G. K., Pickett, H. M., Santee, M. L., Wu,D. L., Boyles, M. A., Burke, J. R., Lay, R. R., Loo, M. S., Livesey,N. J., Lungu, T. A., Manney, G. L., Nakamura, L. L., Perun, V. S.,Ridenoure, B. P., Shippony, Z., Siegel, P. H., Thurstans, R. P., Har-wood, R. S., Pumphrey, H. C. and Filipiak, M. J., 1999: The UARSand EOS Microwave Limb Sounder Experiments. J. Atmos. Sci.,56, 194–218.

Weng, F., 1992: A multi-layer discrete-ordinate method for vector ra-diative transfer in a vertically-inhomogeneous, emitting and scatter-ing atmosphere – I. Theory. J. Quant. Spectrosc. Radiat. Transfer,47, 19–33.