ALIS: An efficient method to compute high spectral resolution polarized solar

radiances using the Monte Carlo approach

Claudia Emde, Robert Buras, Bernhard Mayer

Meteorological Institute, Ludwig-Maximilians-University, Theresienstr. 37, D-80333 Munich, Germany

Abstract
arXiv:1901.01842v1 [astro-ph.EP] 7 Jan 2019

An efficient method to compute accurate polarized solar radiance spectra using the (three-dimensional) Monte Carlo
model MYSTIC has been developed. Such high resolution spectra are measured by various satellite instruments for
remote sensing of atmospheric trace gases. ALIS (Absorption Lines Importance Sampling) allows the calculation
of spectra by tracing photons at only one wavelength: In order to take into account the spectral dependence of the
absorption coefficient a spectral absorption weight is calculated for each photon path. At each scattering event the
local estimate method is combined with an importance sampling method to take into account the spectral dependence
of the scattering coefficient. Since each wavelength grid point is computed based on the same set of random photon
paths, the statistical error is the almost same for all wavelengths and hence the simulated spectrum is not noisy. The
statistical error mainly results in a small relative deviation which is independent of wavelength and can be neglected
for those remote sensing applications where differential absorption features are of interest.
Two example applications are presented: The simulation of shortwave-infrared polarized spectra as measured by
GOSAT from which CO2 is retrieved, and the simulation of the differential optical thickness in the visible spectral
range which is derived from SCIAMACHY measurements to retrieve NO2 . The computational speed of ALIS (for
one- or three-dimensional atmospheres) is of the order of or even faster than that of one-dimensional discrete ordinate
methods, in particular when polarization is considered.
Keywords: radiative transfer, Monte Carlo method, polarization, trace gas remote sensing, high spectral resolution

1. Introduction the individual absorption lines can not be distinguished.

Nevertheless each molecule type has its specific absorp-
Monitoring of atmospheric trace gases is important tion features so that the measured spectra include infor-
to understand atmospheric composition and global cli- mation about the various trace gas concentrations. In
mate change. In particular climate models require in- the near-infrared spectral range there are mainly vibra-
formation about the concentration and global distribu- tional transitions; here individual lines can be identified
tion of trace gases like e.g. H2 O, CO2 , O3 , or CH4 . The and used for trace gas measurements.
trace gases can be observed by measuring solar radia-
tion which is scattered and absorbed by the molecules. Examples for currently operating satellite instru-
Several instruments have been developed: satellite in- ments that measure high resolution radiance spectra of
struments provide global observations, local measure- scattered solar radiation are SCIAMACHY on the EN-
ments can be taken from the ground, from air-plane, VISAT satellite (Gottwald et al. [1]), OMI on AURA
or from a balloon. Most instruments designed for trace (Levelt et al. [2]), GOME-2 on METOP (Callies et al.
gas concentrations observations measure radiance spec- [3]) and TANSO-FTS on GOSAT (Kuze et al. [4]).
tra with high spectral resolution. In the UV-Vis spectral SCIAMACHY and TANSO-FTS have the advantage of
range, absorption of radiation is due to molecular transi- measuring not only the radiance but also the polariza-
tions; at the same time vibrational and rotational transi- tion state of the radiation. While extraterrestrial solar ra-
tions can take place, which results in band spectra where diation is unpolarized, the radiance arriving at the satel-
lite is polarized due to scattering by molecules, aerosols,
or clouds and due to surface reflection. The polarization
Email address: [email protected] (Claudia Emde) information may therefore be used to reduce the uncer-
Preprint submitted to Journal of Quantitative Spectroscopy & Radiative Transfer January 8, 2019
tainties in trace gas retrievals introduced by aerosols, validated ALIS by comparison to the well-known and
clouds and surface reflection. well-tested DISORT program, originally developed and
The retrieval of trace gas concentrations from radi- implemented by Stamnes et al. [7] in FORTRAN77. We
ance spectra requires a so called forward model, which use a new version of the code implemented in C (Buras
can simulate such measurements for given realistic at- et al. [17]) with increased efficiency and numerical ac-
mospheric conditions. For the often used optimal esti- curacy.
mate retrieval method (Rodgers [5]) it is important that
the forward model is fast because it has to be run sev- 2. Methodology
eral times until iteratively the atmospheric composition
is found which best matches the measured spectra. The new method Absorption Lines Importance Sam-
A commonly used method to simulate solar radia- pling (ALIS), which allows fast calculations of spec-
tive transfer is the discrete ordinate method which was tra in high spectral resolution, has been implemented
first described by Chandrasekhar [6] and which has been into the radiative transfer model MYSTIC (Monte Carlo
implemented for instance into the freely available well code for the phYsically correct Tracing of photons In
known software DISORT (Stamnes et al. [7]). The DIS- Cloudy atmospheres; Mayer [18]). MYSTIC is oper-
ORT code however has the limitations that it assumes ated as one of several solvers of the libRadtran radia-
a plane-parallel atmosphere (i.e. horizontal inhomo- tive transfer package (Mayer and Kylling [19]). The
geneities can not be taken into account) and that it ne- common model geometry of MYSTIC is a 3D plane-
glects polarization. Polarization has been included in parallel atmosphere to simulate radiances or irradiances
the VDISORT code (Schulz et al. [8]). The SCIATRAN in inhomogeneous cloudy conditions. The model can
code (Rozanov et al. [9]) is also based on the discrete also be operated in 1D spherical model geometry (Emde
ordinate method. It can optionally take into account and Mayer [20]) which makes it suitable also for limb
spherical geometry as well as polarization (Rozanov and sounding applications. Recently MYSTIC has been ex-
Kokhanovsky [10]). tended to handle polarization due to scattering by ran-
Another method for the simulation of solar radia- domly oriented particles, i.e. clouds, aerosols, and
tive transfer is the Monte Carlo method (Marchuk et al. molecules (Emde et al. [21]), and to handle topography
[11], Marshak and Davis [12]), which is usually much (Mayer et al. [22]). Several variance reduction tech-
slower than the discrete ordinate method. For this rea- niques were also introduced to MYSTIC in order to
son Monte Carlo methods have mostly been used for speed up the computation of radiances in the presence
simulations including inhomogeneous clouds (e.g. Zin- of clouds and aerosols (Buras and Mayer [23]).
ner et al. [13]) for which the plane-parallel approxi-
mation can not be applied. We have developed a new 2.1. Monte Carlo method for solar atmospheric radia-
Monte Carlo method which calculates high spectral res- tive transfer
olution radiance spectra very efficiently. The algorithm, This section briefly describes the implementation of
named ALIS (Absorption Lines Importance Sampling), solar radiative transfer in MYSTIC which is explained
does not require any approximations, in particular it can in detail in Mayer [18]. We describe only those details
easily take into account polarization and horizontal in- which are required to understand the following sections
homogeneity. We show that the computational time of about the ALIS method.
ALIS for high resolution radiance spectra is comparable In the forward tracing mode “photons”1 are traced on
to or even faster than the discrete ordinate approach, es- their way through the atmosphere. The photons are gen-
pecially if polarization is included. This means that the erated at the top of the atmosphere where their initial
algorithm is sufficiently fast to be used as forward model direction is given by the solar zenith angle and the solar
for trace gas retrieval algorithms. The basis of the ALIS azimuth angle.
method is that all wavelengths are calculated at the same Absorption and scattering are treated separately: Ab-
time based on the same random numbers. This method sorption is considered by a photon weight according to
which is sometimes called “method of dependent sam- Lambert-Beer’s law:
Z s
pling” (Marchuk et al. [11]) has been used for various
!
applications, e.g. to calculate mean radiation fluxes in wabs = exp − βabs (s0 )ds0 (1)
0
the near-IR spectral range (Titov et al. [14]), to compute
multiple-scattering of polarized radiation in circumstel- 1 We use the term ”photon” to represent an imaginary discrete
lar dust shells (Voshchinnikov and Karjukin [15]), or to amount of electromagnetic energy transported in a certain direction.
calculate Jacobians (Deutschmann et al. [16]). We have It is not related to the QED photon [24].

2
Here ds0 is a path element of the photon path and βabs = series which is a solution of the integral form of the ra-
i=1 βabs,i is the total absorption coefficient which is the
PN
diative transfer equation (see e.g. Marshak and Davis
sum of the N individual absorption coefficients βabs,i [12]).
of molecules, aerosols, water droplets, and ice crystals. Additional weights are required to take into account
The integration is performed over the full photon path. polarization (Emde et al. [21]) and variance reduction
The free path of a photon until a scattering interac- techniques (Buras and Mayer [23]).
tion takes place is sampled according to the probability After tracing N p photons the radiance is given by the
density function (PDF): average contribution of all photons:
Z s ! PN p
Ii
P(s) = βsca (s) exp − βsca (s )ds
0 0
(2) I= 1 (6)
0 Np

Here βsca = i=1

PN
βsca,i is the total scattering coefficient The methods described above are implemented for
of N interacting particle and molecule types. monochromatic radiative transfer. If one wants to calcu-
We use a random number r ∈ [0, 1] to decide which late a radiance spectrum using these methods one has to
interaction takes place. If there are N types of particles calculate all spectral points subsequently. Here usually
and molecules at the place of scattering, the photon in- a very high accuracy is required in order to distinguish
teracts with the jth type if the random number fulfills the spectral features from statistical noise which means that
following condition: such calculations are computationally expensive.

P j−1 Pj 2.2. Calculation of high spectral resolution clear-sky

βsca,i βsca,i
i=1
<r≤ i=1
(3) radiance spectra
βsca βsca
In the following an efficient method how to com-
At each scattering event the “local estimate” weight pute high spectral resolution radiance spectra will be
is calculated which corresponds to the probability that described and demonstrated on the example of the spec-
the photon is scattered into the direction of the detector tral region from 765–768 nm in the O2 A absorption
and reaches it without further interactions: band where we calculate the spectrum with a resolu-
 R  tion of 0.003 nm. The line-by-line gas absorption co-
P11 (θ p ) exp − (βabs + βsca )ds0 efficients have been computed using the ARTS model
wle,is = (4) (Buehler et al. [25], Eriksson et al. [26]) for the standard
cos(θd )
mid-latitude summer atmosphere (Anderson et al. [27]).
Here θ p is the angle between photon direction (before Fig. 1 shows the vertically integrated optical thickness
scattering) and the radiance direction. The phase func- of molecular absorption τabs,v (top) and scattering τsca,v
tion P11 (first element of the scattering matrix) gives the (bottom). Whereas the scattering optical thickness for
probability that the photon is scattered into the direc- the cloudless, aerosol-free atmosphere is rather small
tion of the detector, “is” denotes the scattering order. In and almost constant, it varies only from about 0.0239
order to calculate the probability that the photon actu- to 0.0243, the absorption optical thickness varies over
ally reaches theRdetector the Lambert-Beer term for ex- five orders of magnitude, from about 10−3 to 102 (note
the logarithmic scale).

tinction exp − (βabs + βsca )ds0 needs to be included.
Finally we have to divide by the zenith angle of the de- As mentioned in the previous section, absorption
tector direction θd to account for the slant area in the is considered separately by calculating the absorption
definition of the radiance. The contribution of the pho- weight wabs (Eq. 1). In order to calculate a radiance
ton to the radiance measured at the detector is then given spectrum taking into account the spectral variation of
as the absorption coefficient βabs we can easily calculate
XNs the absorption weight for each wavelength and get a
Ii = wabs,is wle,is (5) spectrally dependent absorption weight vector:
is=1 Z s !
Here i is the index of the photon, N s is the number of wabs (λ) = exp − βabs (λ, s )ds
0 0
(7)
0
scattering events along the photon path, and wabs,is is
the absorption weight (Eq. 1) evaluated at the scattering Here λ denotes the wavelength of the radiation. In prac-
order is. One can show formally that the sum over the tice the integral corresponding
Rs to the absorption opti-
local estimate weights corresponds to a von Neumann cal thickness τabs = 0 βabs ds0 is calculated step by step
3
 1e+03 0.004

1e+02
0.003

transmittance
1e+01
τabs,v

1e+00 0.002

1e-01
0.001

1e-02

0.000
1e-03 765.0 765.5 766.0 766.5 767.0 767.5 768.0
765.0 765.5 766.0 766.5 767.0 767.5 768.0
1.0
0.0244
0.5

0.0243 0.0

relative difference [%] -0.5

0.0242
-1.0
τsca,v

0.0241 -1.5

-2.0
0.0240
-2.5
0.0239
-3.0

0.0238 -3.5
765.0 765.5 766.0 766.5 767.0 767.5 768.0
wavelength [nm]
765.0 765.5 766.0 766.5 767.0 767.5 768.0
wavelength [nm]
Figure 2: Radiance spectra calculated with MYSTIC in comparison
to DISORT calculations. The top panel shows the transmittance (ra-
Figure 1: Integrated vertical optical thickness of molecular absorption
diance normalized to extraterrestrial irradiance) spectra of two inde-
(top) and molecular scattering (bottom).
pendent MYSTIC runs (grey lines) and the DISORT result (black line)
and the bottom panel shows the relative differences between the MYS-
TIC runs and DISORT in percent.
while the photon moves through the layers/boxes of the
discretized model atmosphere (see Mayer [18]):
viewing angle of 60◦ . We did not include any sensor
X
τabs (λ) = βabs (λ, p)∆s p (8)
p response function. The upper panel shows the transmit-
tance spectra (radiance divided by extraterrestrial irra-
Here the p denotes the step index along the photon path, diance) and the lower panel shows the relative differ-
and ∆s p is the pathlength of step p. We also include ence to the DISORT solver operated with 32 streams.
the spectrally dependent absorption coefficient βabs (λ) The MYSTIC calculations with 106 photons took 13 s
in the local estimate weight wle,is (λ). Thus we only need on a single processor with 2 GHz CPU (all computa-
to trace the photons for one wavelength, calculate the tional times in the following refer to this machine), the
spectral absorption weights and get the full radiances DISORT calculation took 25 s. The relative difference
spectrum with high spectral resolution. For each photon between MYSTIC and DISORT is less than about 2%
we get (compare Eq. 5): with some exceptions where the transmittance is almost
zero. The spectral features in the MYSTIC calcula-
Ns
X tions are well resolved. The two MYSTIC runs used
Ii (λ) = wabs,is (λ)wle,is (λ) (9)
exactly the same setup but the results show a deviation
is=1
between each other and with respect to the DISORT re-
Fig. 2 shows two spectral calculations using this sult. This deviation is due to the statistical error of the
method. Here we assumed that the sun is in the zenith Monte Carlo calculation, with 106 photons the standard
and the sensor is on the ground and measures with a deviation is 0.66%. Hence the deviation can be reduced
4
by running more photons. Since the same photon paths Also, as long as we neglect the wavelength depen-
are used to compute all wavelengths the deviation is the dence of the Rayleigh depolarization factor (see e.g.
same at all spectral data points and the spectra are not Hansen and Travis [29]) the Rayleigh phase function is
noisy. However the deviation shows a spectral depen- wavelength-independent and wsca2,is (λ) = 1.
dence which is not a statistical error but can be attributed The final result for the contribution of a photon in-
to the spectral dependence of Rayleigh scattering which cluding the spectral dependence of absorption and scat-
has been neglected so far. In the calculation βsca was set tering to the spectral radiance is:
to a constant value corresponding to βsca at 766.5 nm.
In the next section we will describe how to include the
Ii (λ) = (14)
spectral dependence of the scattering coefficient.
Ns
X
 is

 Y 
2.3. Importance sampling for molecular scattering wabs,is (λ)  wsca1,is0 (λ, s0 )wsca2,is0 (λ, s0 ) wle,is (λ)
Eq. 2 is the PDF which we use for sampling the free is=1 is0 =1

pathlength of the photons, where the scattering coeffi-

cient βsca now becomes wavelength dependent. We want
Now we calculate again the spectrum in the O2 A-
to use this PDF for sampling the pathlength for all wave-
band from 765–768 nm with λc =765 nm and compare
lengths. In order to ensure that the results are correct for
the result with an accurate DISORT calculation using
all wavelength we introduce a correction weight (impor-
64 streams. The top left panel of Fig. 3 shows the rel-
tance sampling method, see e.g. Ripley [28]):
ative difference of a MYSTIC run with 106 photons,
wsca1 (λ, s) (10) which takes 14.6 s, to DISORT. We see that there is still
 Rs  a relative deviation of about 0.4% which is due to the
βsca (λ, s) exp − 0 βsca (λ, s )ds
0 0
=  Rs  statistical error of the Monte Carlo calculation, but the
βsca (λc , s) exp − 0 βsca (λc , s0 )ds0 spectral dependence of the deviation is now removed
because we have corrected the spectral dependence of
βsca (λ, s)
Z s !
= exp − (βsca (λ, s0 ) − βsca (λc , s0 ))ds0 the scattering coefficient. In order to check whether
βsca (λc , s) 0 the method is correctly implemented without any bias
Here λc (c for “computational”) is the wavelength cor- (apart from the statistical error) we performed a MYS-
responding to the scattering coefficient that is used to TIC run with 109 photons. The result is shown in the
sample the photon free path. As in the previous section lower left panel. The spectrally independent deviation
we write the second part of this expression as a sum over has almost vanished (<0.01%) and the relative differ-
the model layers/boxes: ence between MYSTIC and DISORT is below 0.05%.
βsca (λ, s) For comparison we show in the right panels of the figure
wsca1 (λ, s) = (11) DISORT runs with 16 and 32 streams respectively com-
βsca (λc , s)
  pared to the DISORT run with 64 streams. The differ-
ence between DISORT (16 streams) and DISORT (64
 X 
× exp − (βsca (λ, p) − βsca (λc , p))∆s p 
p
streams) is actually larger than the difference between
MYSTIC (109 photons) and DISORT (64 streams).
The probability that the photon is scattered into a di-
It should be noted that this Monte Carlo method does
rection with scattering angle θ p is given by the phase
only work well as long as the scattering coefficient does
function P11 (λ, θ p ). So we need another weight to cor-
not vary too much within the simulated wavelength re-
rect for the spectral dependence of the phase function
gion. Else the scattering weights can obtain values very
which can again easily be derived using the importance
far from 1, resulting in large statistical noise and slow
sampling method:
convergence.
P11 (λ, θ p , s)
wsca2,is (λ, s) = (12)
P11 (λc , θ p , s)
Here s is the location at which the photon is scattered. 2.4. Calculation of high resolution spectra including
Note that in the case where we have only molecules as aerosol and cloud scattering
interacting particles and neglect depolarization P11 is
the Rayleigh phase function
3 It is straightforward to apply the method to an atmo-
P11 (θ p ) = (1 + cos2 θ p ) (13) sphere including clouds and/or aerosols. We just need
4
5
 MYSTIC with 106 photons DISORT with 16 streams
relative difference [%] 1.0 0.10

0.5 0.05

0.0 0.00

−0.5 −0.05

−1.0 −0.10
765.0 765.5 766.0 766.5 767.0 767.5 768.0 765.0 765.5 766.0 766.5 767.0 767.5 768.0

MYSTIC with 10 photons

9
DISORT with 32 streams
0.10 0.10
relative difference [%]

0.05 0.05

0.00 0.00

−0.05 −0.05

−0.10 −0.10
765.0 765.5 766.0 766.5 767.0 767.5 768.0 765.0 765.5 766.0 766.5 767.0 767.5 768.0
wavelength [nm] wavelength [nm]

Figure 3: Relative differences of various model setups with respect to a DISORT calculation with 64 streams. The left panels show MYSTIC
calculations with 106 and 109 photons respectively. The right panels show DISORT calculations with 16 and 32 streams respectively.

to use the total absorption and scattering coefficients 3. Applications

3.1. Simulation of polarized near infrared spectra in
N
X cloudless conditions
βabs (λ) = βabs,i (λ) (15) The Greenhouse Gases Observing Satellite (GOSAT)
i=1 determines the concentrations of carbon dioxide and
N
X methane globally from space. The spacecraft was
βsca (λ) = βsca,i (λ) (16) launched on January 23, 2009, and has been op-
i=1 erating properly since then. Information about the
project can be found on the web-page http://www.
gosat.nies.go.jp. GOSAT carries the Thermal and
and the average phase function given by
Near Infrared Sensor for Carbon Observation Fourier-
Transform Spectrometer (TANSO-FTS) (Kuze et al.
PN
βsca,i (λ)P11,i (λ) [4]) which measures in 4 spectral bands (band 1: 0.758–
P11 (λ) = i=1
(17) 0.775 µm, band 2: 1.56–1.72 µm, band 3: 1.92–
βsca (λ)
2.08 µm, band 4: 5.56–14.3 µm). The spectral reso-
lution in all bands is 0.2 cm−1 . For the visible spec-
Here N is the number of interacting particles/molecules. tral range (bands 1–3) polarized observations are per-
These quantities can be used to compute the wavelength formed. In order to analyze this data a fast polar-
dependent weights wabs (λ) (Eq. 7), wsca1 (λ) (Eq. 11) and ized radiative transfer code is required. The Monte
wsca2 (λ) (Eq. 12). In MYSTIC we so far consider only Carlo approach which is described in this study is an
the spectral dependence of molecular scattering because alternative to commonly used discrete ordinate or dou-
the spectral dependence of cloud and aerosol scattering bling and adding codes. The approach is fully com-
can safely be neglected in narrow wavelength intervals. patible to the implementation of polarization in MYS-
6
 I
TIC as described in Emde et al. [21] and validated 0.06
in Kokhanovsky et al. [30], because the weight vector
which is calculated to take into account the polarization 0.05

state of the photon does not interfere with the spectral

0.04
weights. An advantage of the Monte Carlo approach is

reflectance
of course that it is easy to take into account horizontal 0.03
inhomogeneities of clouds, aerosols, and molecules.
In the following we show an example simulation 0.02
where we selected a spectral window of band 3 from
0.01
4815–4819 cm−1 (corresponding to ≈2.075–2.077 µm)
in the near infrared. The radiance simulation is per-
0.00
formed with a spectral resolution of 0.01 cm−1 . The at- 4815 4816 4817 4818 4819

mospheric profiles (pressure, temperature, trace gases) Q

0.000
correspond to the standard mid-latitude summer atmo-
sphere as defined by Anderson et al. [27]. The molec-
-0.005
ular absorption coefficients have been computed using
the ARTS line-by-line model. We assume a thin mar- reflectance
-0.010
itime aerosol layer with an optical thickness of 0.05
at 2 µm. We took the refractive indices and the size
distribution data from Hess et al. [31] (maritime clean -0.015

aerosol mixture) and performed Mie calculations to ob-

tain the aerosol optical properties including the phase -0.020

matrix. We assume an underlying water surface which

is modeled using the reflectance matrix as defined in -0.025
4815 4816 4817 4818 4819
Mishchenko and Travis [32]. The reflectance matrix wavenumber [cm−1 ]
is based on the Fresnel equations, on Cox and Munk
[33, 34] to describe the wind-speed dependent slope of Figure 4: Simulated GOSAT spectra over ocean. The assumed wind
speed is 5 m/s. The viewing angle of the FTS is 30◦ . The line-styles
the waves, and on Tsang et al. [35] to account for shad- correspond to different solar zenith angles, the solid line corresponds
owing effects. The wind speed was taken to be 5m/s. to 30◦ , i.e. the sun glint is observed, the dashed and the dash-dotted
The viewing angle of the FTS is 30◦ and simulations lines correspond to 20◦ and 60◦ , respectively. All simulations are in
have been performed for the principal plane assuming the principal plane. The upper panel shows the normalized reflected
intensity I and the lower panel shows the polarization difference Q.
different solar zenith angles θ0 . The full Stokes vector
has been computed for all setups.
Fig. 4 shows the simulated GOSAT spectra. The solid measured spectra are converted to so-called differential
lines correspond to θ0 =30◦ , in this case the FTS ob- optical thicknesses defined as
serves the center of the sun glint, therefore this spec-
trum shows the highest reflectance. The dashed lines are D(λ) = ln(ITOA (λ)) − P3 (λ) (18)
for θ0 =20◦ , still in the sun-glint region and the dashed-
where ITOA (λ) is the reflectance at the top of the at-
dotted lines are for θ0 =60◦ which is not influenced much
mosphere. P3 (λ) is a third degree least square polyno-
by the sun glint. The computation time for each po-
mial fit of the logarithm of ITOA (λ) with respect to the
larized spectrum using 106 photons was 2 minutes and
wavelength, which removes the slowly varying part of
25 seconds, the standard deviation (approximately the
the spectrum. The conversion of ITOA (λ) into D(λ) im-
same for each Stokes vector component) for θ0 =20◦ and
proves the contrast of the NO2 absorption line depths
θ0 =30◦ is 0.03%, for θ0 =60◦ it is 0.16%.
and thereby the accuracy of the retrieval. The retrieval
algorithm minimizes the function
3.2. Simulations of differential optical thickness in bro-  
ken cloud conditions F(λ, VNO2 ,ret ) = D(λ, VNO2 ,true ) − D(λ, VNO2 ,ret ) (19)

Retrievals of the tropospheric NO2 column from where VNO2 ,true and VNO2 ,ret are the true and the retrieved
SCIAMACHY data are based on the differential optical tropospheric NO2 columns, respectively.
absorption spectroscopy (DOAS) method (Richter and Our new method allows efficient computations of
Burrows [36], Richter et al. [37]). For this method the D(λ). As an example Fig. 5 shows spectra for three
7
 0.003
ference is less than 0.1% and shows the typical Monte
0.002 Carlo noise.
The bottom left panel shows the differential optical
0.001
differential optical depth

thicknesses derived from the simulations. Here the sta-

tistical noise of the MYSTIC calculation without ALIS
0.000
(grey line) is clearly visible. All ALIS simulations re-
-0.001 sult in a smooth differential optical thickness because
all wavelengths are calculated based on the same pho-
-0.002 ton paths. This yields a relative deviation in the re-
flectance which is independent of wavelengths and can
-0.003
be removed completely by subtracting the fitted poly-
-0.004 nomial. The bottom middle panel shows the absolute
difference between the differential optical thicknesses
-0.005 derived from the ALIS simulations and the one derived
400 410 420 430 440 450 460 470
wavelength [nm] from the DISORT simulation. Even for the simulation
with only 103 photons the differential optical thickness
Figure 5: Differential optical thickness calculated for three different
NO2 profiles, corresponding to low (- -), medium (· · · ) and highly pol-
is quite accurate, the difference w.r.t. DISORT is of the
luted (—) conditions. The spectra have been computed using MYS- order of a few per-cent. Using 105 photons or more
TIC with 105 photons. yields very accurate differential optical thicknesses, the
difference is here well below 1%. The bottom right
panel shows the difference between the MYSTIC calcu-
different NO2 profiles, corresponding to low, medium lation without ALIS and the DISORT calculation. The
and highly polluted conditions. The Lambertian surface difference is of the same order of magnitude as the dif-
albedo was set to 0.1 and the solar zenith angle to 32◦ . ferential optical thickness itself and hence the accuracy
NO2 and O3 profiles are the same as used in the study of this simulation would not be sufficient for NO2 re-
by Vidot et al. [38]. The NO2 absorption cross sections trievals.
have been taken from Burrows et al. [39], ozone absorp- Fig. 6 clearly demonstrates that a common Monte
tion was also included in the simulations using the cross Carlo approach which calculates wavelength by wave-
sections by Molina and Molina [40]. The spectral reso- length sequentially is extremely inefficient for simula-
lution of the simulation is 0.1 nm. tions of the differential optical thickness because each
Fig. 6 shows calculations for the NO2 profile corre- wavelength has an independent statistical error which is
sponding to medium polluted conditions. The top left larger than the absorption features unless a huge number
panel shows the reflectance, where the grey line corre- of photons is used. In order to obtain a result with an
sponds to a MYSTIC calculation without using ALIS, accuracy comparable to ALIS with 103 photons (0.9 s
i.e. all wavelengths are simulated subsequently using computation time), MYSTIC without ALIS would re-
107 photons for each wavelength. The standard devia- quire at least 109 photons per wavelength which would
tion for each wavelength is about 0.03%. This calcula- take about 138 days computation time.
tion took 33 h 12 min on one CPU. The black line shows The impact of cirrus clouds on tropospheric NO2 re-
the calculations using ALIS with different numbers of trievals has been investigated in a sensitivity study by
photons. The calculation using 103 photons took 0.9 s, Vidot et al. [38]. They take into account the sub-pixel
the one with 105 photons took 38 s, and the one with inhomogeneity by a simple independent pixel approach.
107 photons took 63 min 53 s. Obviously the absorption Using our Monte Carlo approach we can take into ac-
features of NO2 are barely visible in the reflectance plot. count full 3-dimensional radiative transfer, e.g. the in-
There is a deviation between the simulations which is teraction between the cloudy and the clear-sky part of
due to the statistical error of the simulation using ALIS the domain. Fig. 7 shows an example where the setup
with 103 or 105 photons. The top middle panel shows is similar to the study by Vidot et al. [38]. We have
the relative differences of the ALIS simulations w.r.t. taken a very simple 3D cloud field, the cirrus clouds
DISORT, which requires 30 s computation time with were modeled as 1×1×1 km3 cubes and arranged as a
32 streams. Obviously the relative difference decreases chess-board, hence the cloud fraction c is 0.5. The sur-
with increasing number of photons. The top right panel face albedo is 0.05 and the solar zenith angle is 30◦ .
shows the relative difference between the MYSTIC cal- The optical thickness of the clouds is 3, the geometrical
culation without ALIS and DISORT. The relative dif- thickness is 1 km, the cloud base height is 10 km and
8
 0.054
MYSTIC lbl ALIS 103 ph. MYSTIC lbl
0.052 4.0 0.4
ALIS 103 ph. ALIS 105 ph.

relative difference [%]

relative difference [%]

0.050 ALIS 105 ph. ALIS 107 ph.
2.0 0.2
ALIS 107 ph.
reflectance

0.048

0.046 0.0 0.0

0.044
-2.0 -0.2
0.042

0.040 -4.0 -0.4

0.038
400 410 420 430 440 450 460 470 400 410 420 430 440 450 460 470 400 410 420 430 440 450 460 470

2.0e-03 1.0e-03
1.5e-03 4.0e-05
differential optical depth

1.0e-03
absolute difference

absolute difference
5.0e-04
2.0e-05
5.0e-04
0.0e+00
0.0e+00 0.0e+00
-5.0e-04
-1.0e-03
-2.0e-05
-1.5e-03 -5.0e-04

-2.0e-03 -4.0e-05
-2.5e-03 -1.0e-03
400 410 420 430 440 450 460 470 400 410 420 430 440 450 460 470 400 410 420 430 440 450 460 470
wavelength [nm] wavelength [nm] wavelength [nm]

Figure 6: Simulations for NO2 profile corresponding to medium polluted conditions. The top left panel shows the reflectance and the bottom left
panel the differential optical thickness. The black lines correspond to Monte Carlo simulations using the ALIS method with different number of
photons (103 , 105 , and 107 ), the grey line shows a Monte Carlo simulation without ALIS (all wavelengths are calculated independently). The
middle panels show differences w.r.t. DISORT for the ALIS simulations and the right plots show differences w.r.t. DISORT for the MYSTIC
calculation without using ALIS.

the effective radius is 30 µm where the parameterization calculation, because photons which are scattered out of
by Baum et al. [41, 42] was used for the cirrus optical the clouds on the sides have a higher probability of be-
properties. The solar zenith angle is 30◦ and the wave- ing transmitted to the surface.
length range is 400 nm to 470 nm. We performed 3D The bottom panel of Fig. 7 shows the differential op-
calculation and also used the independent pixel approx- tical thickness. The difference between IPA and 3D is in
imation (IPA) for comparison. All simulations shown in this case about 10% which will cause an error of some
Fig. 7 were calculated using MYSTIC with ALIS. The per cent in the tropospheric NO2 retrievals. Note that
reflectance for the IPA simulation is calculated as the this calculation is only an example to demonstrate the
sum of the reflectance of the clear-sky part Rclear and the new algorithm to calculate high spectral resolution spec-
reflectance of the cloudy part Rcloud weighted with the tra using the Monte Carlo method. With different setups
cloud fraction: the error on the retrieval can be completely different.

R = cRcloud + (1 − c)Rclear (20)

4. Conclusions
In order to speed up the calculations in the presence
of clouds, the variance reduction technique VROOM We have developed the new method ALIS (absorp-
(Buras and Mayer [23]) was used. Using VROOM the tion lines importance sampling) that allows to compute
simulation 105 photons are sufficient to obtain an accu- polarized radiances in high spectral resolution using the
rate result with a standard deviation of approximately Monte Carlo method in a very efficient way. We sample
0.5%. The 3D calculation using these settings took random photon paths at one wavelength. For these ran-
1 min 34 s, an IPA calculation using DISORT with 32 dom paths we calculate a spectral absorption weight us-
streams takes 58 s. ing the wavelength dependent absorption coefficients of
Fig. 7 shows a part of the spectrum where we have the model boxes. In order to correct for the wavelength
pronounced features in the differential optical depth. In dependence of Rayleigh scattering an importance sam-
the top panel one can see that for all wavelengths the re- pling method is applied. If necessary the same method
flectance in the 3D calculation is smaller than in the IPA can be applied to correct for the spectral dependence
9
 0.09
ferential absorption features are of interest, this devia-
0.08 tion does not matter.
Two example applications are shown. First the simu-
0.07 lation of polarized near-infrared spectra over an ocean
reflectance

surface as measured by e.g. GOSAT. Here we simu-

0.06
lated the Stokes vector with a standard deviation smaller
0.05 than 0.05% for 400 spectral points in 2 minutes 25 sec-
onds on a single PC. For a standard deviation of 0.5%
0.04 the calculation would be 100 times faster. These short
computation times show that the algorithm has the po-
0.03
tential to be used as a forward model for trace gas re-
0.02 trievals from polarized radiance measurements, in par-
430 435 440 445 450
ticular since commonly used discrete ordinate methods
become much slower when they include polarization.
0.003
The second example is the simulation of the differ-
0.002 ential optical thickness from 400 nm to 470 nm which
differential optical depth

is used to retrieve NO2 from e.g. SCIAMACHY. Here

0.001 the computation time for accurate (scalar) simulations
0.000
was comparable to DISORT. We performed this calcu-
lation also for an inhomogeneous cloud scene where
-0.001 cirrus clouds are approximated by simple cubes. We
compared the result of the 3D simulation with an in-
-0.002 dependent pixel calculation and found a difference of
-0.003
about 10 % in the differential optical thickness for this
example. The calculations show that ALIS is suitable to
-0.004 study effects of horizontal inhomogeneity on trace gas
430 435 440 445 450
wavelength [nm] retrievals in presence of cirrus clouds.

Figure 7: Impact of cirrus clouds on a spectrum used for NO2 re-

trievals. The cirrus clouds are modeled as 1×1×1 km3 cubes, the
Acknowledgement
cloud fraction in the model domain is 0.5. Black lines correspond
to the independent pixel (IPA) calculation and grey lines to the 3D We thank Timothy E. Dowling for translating the
calculation. Upper panel: The thick lines show the reflectance R of DISORT code from FORTRAN77 to C which resulted
the full domain. The dashed line shows the the clear-sky pixels only in great improvements regarding numerical accuracy
(Rclear ) and the thin solid line shows the cloudy pixels (Rcloud ) for the and computation time. Furthermore we thank Jerôme
IPA calculation. The lower panel shows corresponding differential
optical thicknesses, where the lines styles are defined as above. Vidot for providing NO2 profiles. This work was done
within the project RESINC2 funded by the “Deutsche
Forschungsgemeinschaft” (DFG).
of cloud and aerosol scattering. The method allows
us to calculate radiances for many wavelengths at the
References
same time without significantly increasing the compu-
tational cost. ALIS has been implemented in the MYS- References
TIC model which handles three-dimensional radiative
[1] M. Gottwald, H. Bovensmann, G. Lichtenberg, S. Noel, A. von
transfer in cloudy atmospheres including polarization Bargen, S. Slijkhuis, A. Piters, R. Hoogeveen, C. von Savi-
and complex topography. gny, M. Buchwitz, A. Kokhanovsky, A. Richter, A. Rozanov,
The new algorithm ALIS has been validated by com- T. Holzer-Popp, K. Bramstedt, J.-C. Lambert, J. Skupin, F. Wit-
parison to the well-known and well-tested DISORT trock, H. Schrijver, J. Burrows, SCIAMACHY, Monitoring the
Changing Earth’s Atmosphere, DLR, 2006.
solver. It has been shown that ALIS does not produce [2] P. F. Levelt, G. H. J. van den Oord, M. R. D. adm A. Mlkki,
any bias apart from the statistical noise. Since all wave- H. Visser, J. de Vries, P. Stammes, J. Lundell, H. Saari,
lengths are computed at once, the statistical error is the The Ozone Monitoring Instrument, IEEE Transactions on Geo-
science and Remote Sensing 44 (5) (2006) 1093–1101.
same at all wavelength which results mainly in a relative
[3] J. Callies, E. Corpaccioli, M. Eisinger, A. Hahne, A. Lefebvre,
deviation which is independent of wavelength. How- GOME-2 - Metop’s Second Generation sensor for Operational
ever, for remote sensing applications, where mostly dif- Ozone Monitoring, ESA Bulletin 102.

10
 [4] A. Kuze, H. Suto, M. Nakajima, T. Hamazaki, Thermal and tions from the complex topography of Arizona’s Meteor Crater,
near infrared sensor for carbon observation Fourier-transform Atmos. Chem. Phys. 10 (18) (2010) 8685–8696.
spectrometer on the Greenhouse Gases Observing Satellite for [23] R. Buras, B. Mayer, Efficient unbiased variance reduction tech-
greenhouse gases monitoring, Appl. Opt. 48 (2009) 6716–6733. niques for Monte Carlo simulations of radiative transfer in
[5] C. Rodgers, Inverse Methods for Atmospheric Sounding: The- cloudy atmospheres: The solution, J. Quant. Spectrosc. Radiat.
ory and Practice, World Scientific, 2000. Transfer 112 (3) (2011) 434–447.
[6] S. Chandrasekhar, Radiative transfer, Oxford Univ. Press, UK, [24] M. I. Mishchenko, Gustav Mie and the fundamental concept of
1950. electromagnetic scattering by particles: A perspective, J. Quant.
[7] K. Stamnes, S.-C. Tsay, W. Wiscombe, I. Laszlo, DISORT, Spectrosc. Radiat. Transfer 110 (14-16) (2009) 1210–1222.
a General-Purpose Fortran Program for Discrete-Ordinate- [25] S. A. Buehler, P. Eriksson, T. Kuhn, A. von Engeln, C. Verdes,
Method Radiative Transfer in Scattering and Emitting Layered ARTS, the atmospheric radiative transfer simulator, J. Quant.
Media: Documentation of Methodology, Tech. Rep., Dept. of Spectrosc. Radiat. Transfer 91 (1) (2005) 65–93.
Physics and Engineering Physics, Stevens Institute of Technol- [26] P. Eriksson, S. Buehler, C. Davis, C. Emde, O. Lemke, ARTS,
ogy, Hoboken, NJ 07030, 2000. the atmospheric radiative transfer simulator, version 2, J. Quant.
[8] F. M. Schulz, K. Stamnes, F. Weng, VDISORT: An improved Spectrosc. Radiat. Transfer In Press, Uncorrected Proof.
and generalized discrete ordinate method for polarized (vector) [27] G. Anderson, S. Clough, F. Kneizys, J. Chetwynd, E. Shet-
radiative transfer, J. Quant. Spectrosc. Radiat. Transfer 61 (1) tle, AFGL atmospheric constituent profiles (0-120 km), Tech.
(1999) 105–122. Rep. AFGL-TR-86-0110, Air Force Geophys. Lab., Hanscom
[9] A. Rozanov, V. Rozanov, M. Buchwitz, A. Kokhanovsky, J. P. Air Force Base, Bedford, Mass., 1986.
Burrows, SCIATRAN 2.0 A new radiative transfer model for [28] B. D. Ripley, Stochastic Simulation, Wiley, New York, 2006.
geophysical applications in the 175–2400 nm spectral region, [29] J. E. Hansen, L. D. Travis, Light scattering in planetary atmo-
Adv. Space Res. 36 (2005) 1015–1019. spheres, Space Science Reviews 16 (1974) 527–610.
[10] V. V. Rozanov, A. A. Kokhanovsky, The solution of the vec- [30] A. A. Kokhanovsky, V. P. Budak, C. Cornet, M. Duan, C. Emde,
tor radiative transfer equation using the discrete ordinates tech- I. L. Katsev, D. A. Klyukov, S. V. Korkin, L. C-Labonnote,
nique: Selected applications, Atmospheric Research 79 (2006) B. Mayer, Q. Min, T. Nakajima, Y. Ota, A. S. Prikhach, V. V.
241–265. Rozanov, T. Yokota, E. P. Zege, Benchmark results in vector at-
[11] G. I. Marchuk, G. A. Mikhailov, M. A. Nazaraliev, The Monte mospheric radiative transfer, J. Quant. Spectrosc. Radiat. Trans-
Carlo methods in atmospheric optics, Springer Series in Optical fer 111 (12-13) (2010) 1931–1946.
Sciences, Berlin: Springer, 1980. [31] M. Hess, P. Koepke, I. Schult, Optical Properties of Aerosols
[12] A. Marshak, A. Davis, 3D Radiative Transfer in Cloudy Atmo- and Clouds: The Software Package OPAC, Bulletin of the
spheres, Springer, iSBN-13 978-3-540-23958-1, 2005. American Meteorological Society 79 (5) (1998) 831–844.
[13] T. Zinner, G. Wind, S. Platnick, A. S. Ackerman, Testing re- [32] M. I. Mishchenko, L. D. Travis, Satellite retrieval of aerosol
mote sensing on artificial observations: impact of drizzle and properties over the ocean using polarization as well as intensity
3-D cloud structure on effective radius retrievals, Atmos. Chem. of reflected sunlight, J. Geophys. Res. 102 (1997) 16989–17013.
Phys. 10 (19) (2010) 9535–9549. [33] C. Cox, W. Munk, Measurement of the roughness of the sea sur-
[14] G. Titov, T. B. Zhuravleva, V. E. Zuev, Mean radiation fluxes in face from photographs of the sun’s glitter, Journal of the Optical
the near-IR spectral range: Algorithms for calculation, J. Geo- Society of America 44 (11) (1954) 838–850.
phys. Res. 102 (D2) (1997) 1819–1832. [34] C. Cox, W. Munk, Statistics of the sea surface derived from sun
[15] N. V. Voshchinnikov, V. V. Karjukin, Multiple scattering of po- glitter, Journal of Marine Research 13 (1954) 198–227.
larized radiation in circumstellar dust shells, Astron. Astrophys. [35] L. Tsang, J. A. Kong, R. T. Shin, Theory of Microwave Remote
288 (1994) 883–896. Sensing, John Wiley, New York, 1985.
[16] T. Deutschmann, S. Beirle, U. Frie, M. Grzegorski, C. Kern, [36] A. Richter, J. P. Burrows, Tropospheric NO2 from GOME mea-
L. Kritten, U. Platt, C. Prados-Román, J. Puķite, T. Wagner, surements, Advances in Space Research 29 (11) (2002) 1673 –
B. Werner, K. Pfeilsticker, The Monte Carlo atmospheric radia- 1683.
tive transfer model McArtim: Introduction and validation of Ja- [37] A. Richter, J. P. Burrows, H. N, C. Granier, U. Niemeier, In-
cobians and 3D features, Journal of Quantitative Spectroscopy crease in tropospheric nitrogen dioxide over China observed
and Radiative Transfer 112 (6) (2011) 1119 – 1137. from space, Nature 437 (2005) 129–132.
[17] R. Buras, T. Dowling, C. Emde, New secondary-order inten- [38] J. Vidot, O. Jourdan, A. A. Kokhanosvky, F. Szczap, V. Giraud,
sity correction in DISORT with increased efficiency for for- V. V. Rozanov, Retrieval of tropospheric NO2 columns from
ward scattering, J. Quant. Spectrosc. Radiat. Transfer Submitted satellite measurements in presence of cirrus: A theoretical sen-
2011. sitivity study using SCIATRAN and prospect application for the
[18] B. Mayer, Radiative transfer in the cloudy atmosphere, Euro- A-Train, J. Quant. Spectrosc. Radiat. Transfer 111 (4) (2010)
pean Physical Journal Conferences 1 (2009) 75–99. 586–601.
[19] B. Mayer, A. Kylling, Technical note: The libRadtran software [39] J. P. Burrows, A. Dehn, B. Deters, S. Himmelmann, A. Richter,
package for radiative transfer calclations – description and ex- S. Voigt, J. Orphal, Atmospheric remote–sensing reference data
amples of use, Atmos. Chem. Phys. 5 (2005) 1855–1877. from GOME: Part 1. Temperature–dependent absorption cross
[20] C. Emde, B. Mayer, Simulation of solar radiation during a total sections of NO2 in the 231–794 nm range, J. Quant. Spectrosc.
solar eclipse: A challenge for radiative transfer, Atmos. Chem. Radiat. Transfer 60 (1998) 1025–1031.
Phys. 7 (2007) 2259–2270. [40] L. T. Molina, M. J. Molina, Absolute absorption cross sections
[21] C. Emde, R. Buras, B. Mayer, M. Blumthaler, The impact of of ozone in the 185 to 350 nm wavelength range, J. Geophys.
aerosols on polarized sky radiance: model development, valida- Res. 91 (1986) 14501–14508.
tion, and applications, Atmos. Chem. Phys. 10 (2) (2010) 383– [41] B. Baum, A. Heymsfield, P. Yang, S. Bedka, Bulk scattering
396. models for the remote sensing of ice clouds. Part 1: Micro-
[22] B. Mayer, S. W. Hoch, C. D. Whiteman, Validating the MYS- physical data and models, J. of Applied Meteorology 44 (2005)
TIC three-dimensional radiative transfer model with observa- 1885–1895.

11
[42] B. Baum, P. Yang, A. Heymsfield, S. Platnick, M. King, Y.-X.
Hu, S. Bedka, Bulk scattering models for the remote sensing of
ice clouds. Part 2: Narrowband models, J. of Applied Meteorol-
ogy 44 (2005) 1896–1911.

12