MOCRA: a Monte Carlo code for the simulation of radiative transfer in the atmosphere

Margherita Premuda; Elisa Palazzi; Fabrizio Ravegnani; Daniele Bortoli; Samuele Masieri; Giorgio Giovanelli

doi:10.1364/OE.20.007973

1. Introduction

Propagation of light in the atmosphere can be mathematically described by the Radiative Transfer Equation (RTE), whose solution can be obtained both analytically [1,2] and numerically [3]. The atmosphere is a non homogeneous medium where absorption, scattering and polarization effects combine themselves and make the system extremely difficult for analytical treatment. The Monte Carlo (MC) method is a powerful numerical tool for performing radiative transfer simulations including absorption, polarization and multiple scattering. The method requires the knowledge of the atmospheric physical properties, such as the interaction coefficients and scattering phase functions, obtained by means of different theoretical models according to the type of atmospheric constituents (gases or aerosol particles). Several examples of MC studies of radiative transfer in scattering and absorbing media can be found in the literature [4–6].

The Monte Carlo approach for atmospheric radiative transfer consists in using probabilistic methods to trace the trajectories of individual photons, which are subject to random events. A random number, for example, can be used to decide whether a photon gets scattered or absorbed on its way, and to fix the flight direction in case of scattering. The strict Monte Carlo experiment consists in firing photons out of a radiation source along an initial direction and in following their path through the atmosphere where they are scattered, by which compound and in which direction, till they hit the detector or they are absorbed. The problem in this “forward” procedure is that simulating realistic dimensions and aperture angle of a detector can result in a very low probability for one single photon to get recorded by the modeled receiving apparatus. This would require a very high number of photon histories to be processed to obtain statistically significant results. To overcome this problem, because of linearity of the RTE and reciprocity relationship between the Green’s function and its adjoint [7,8], it is possible to write an adjoint RTE and to perform a “backward” Monte Carlo simulation. Each photon history is traced backward through the atmosphere as being generated by the detector in the line of sight direction, with a photon “weight” proportional to the source intensity, and coming out of the atmosphere towards the sun. Examples of such a simulation scheme can be found in references [9, 10], while in references [4–6] the forward scheme is followed. It is worth pointing out that, in case of atmospheric simulations, a backward Monte Carlo scheme is even simpler than in other contexts, as for instance photon transport in shielding calculations for nuclear reactors. Photons interacting with atmospheric constituents, in fact, do not undergo energy changes, which would require specific and careful treatment [8]. It has been demonstrated that the backward Monte Carlo technique is superior to the forward method in real case applications to atmospheric problems, by greatly reducing the computational demands [11].

A certain number of backward Monte Carlo RTMs have been developed to simulate remote sensing observations of trace gas amounts measured through the Differential Optical Absorption Spectroscopy (DOAS) technique [12–16]. DOAS is an established remote sensing method which identifies and quantifies the trace gases in the atmosphere, taking advantage of their absorption structures in the near UV and visible wavelengths of the solar spectrum (DOAS measuring solar or other natural radiation is called “passive DOAS”). The molecular absorption is analyzed to obtain the concentration of trace gases integrated along the optical path (or Slant Column Density, SCD) between the sun and the receiver. In the case of a vertical looking ground-based sensor (“zenith” configuration), radiative transfer models are used to convert the measured SCD into Vertical Column Density (the concentration of the absorber integrated along the vertical, VCD), via the calculation of the so-called Air Mass Factor (AMF). Perliski and Solomon [17], for instance, developed a multiple scattering backward Monte Carlo RTM to calculate the AMF of their target trace gases from DOAS observations performed with an upward-looking detector. They compared their results with those obtained with a single scattering model [18], and concluded that multiple scattering significantly affects the measurements at twilight. The TRACY (Trace gas RAdiative Transfer Monte Carlo Y(I)implementation) series models [19,20] solve the radiative transfer equation using the backward Monte Carlo method and derive so-called “box air mass factors” which describe the sensitivity of the measurements as a function of atmospheric layer altitude. The TRACY-II version [21], in particular, has been used to retrieve the vertical profiles of the absorbing trace gases from limb-viewing satellite measurements.

More recently, DOAS has been deployed also in “off-axis” configuration [22–24], i.e. ground-based instrument with viewing direction other than vertical one, for a better characterization of tropospheric absorbers. A very recent application of the off-axis configuration is the FRE-DOAS (Flow Rate Emission with DOAS) technique, which was used by Premuda et al. [25] to calculate gaseous pollutant flow rates from moving sources like ships. The AMF concept can be generalized to off-axis configuration measurements and appropriate RTMs have to be used for data interpretation. Palazzi et al. [26] developed the PROMSAR (PROcessing of Multi-Scattered Atmospheric Radiation) code to perform backward Monte Carlo simulations of DOAS observations for vertical and off-axis looking detectors. PROMSAR Box AMFs and radiances were validated through comparison with a series of state-of-the-art UV/visible RTMs [27].

In this paper, the backward Monte Carlo model MOCRA (MOnte Carlo Radiance Analysis) is described. MOCRA is an up-to-date version of the PROMSAR model [26] derived from the PREMAR code [28,29]. MOCRA performs radiative transfer calculations in the UV/Visible spectral range, with the possibility to simulate complex topography, surface albedos and a highly inhomogeneous atmosphere, both vertically and horizontally. The model is suitable for the simulation of different observational set-up, including satellite-limb and nadir viewing geometries, which were not foreseen in the PROMSAR model.

To evaluate the effects on radiative transfer of small variations in some atmospheric constituents, so-called “perturbative” calculations can be performed following photon histories in several perturbed environments at a time. This makes it possible to evaluate a number of effects that could be masked by the variance of the calculations if separate simulations are performed. MOCRA is based on the general scheme of a backward Monte Carlo simulation, taking into account all important physical phenomena except polarization and applying user-selected end of history criteria and variance reducing techniques. It requires the availability of a library data set, being the source of all physical data (vertical spatial distributions of the interaction coefficients, phase functions, refraction indexes etc.) needed for the simulation. For historical reasons related to the development of the previous versions of the MOCRA code, the library source is built running a modified version of the MODTRAN4 code to extract the (macroscopic) interaction coefficients of the atmospheric constituents and their vertical distribution. Extracting the microscopic cross sections for absorption and scattering from the MODTRAN code turned out to be more complicated than for the microscopic interaction coefficient. We plan to use other atmospheric inputs, among which user-defined libraries, than that provided by the MODTRAN code in the future versions of the MOCRA code. This is due to the fact that the library source is historically a modified version of the LOWTRAN-MODTRAN series codes in which interaction coefficients are easier to extract then microscopic cross sections. Making it easier the building of the library for the user, including only microscopic cross sections and vertical distributions of gas and particle densities, will be a future improvement of both PREMAR and MOCRA. When several types of atmosphere have to be taken into account simultaneously in the same calculation, as it occurs when different atmospheric profiles are foreseen for different horizontal regions or in perturbative calculations, the physical data of each region are opportunely put together in an unique input library. Results presented in this paper have been obtained using atmospheric libraries built up from the output of the state-of-the-art MODTRAN4 (MODerate spectral resolution atmospheric TRANsmittance algorithm and computer model) model [30–33], a commercially available RTM developed by the Air Force Research Labs (AFRL) in collaboration with Spectral Sciences, Inc (SSI). MODTRAN is a reference radiation transfer algorithm widely used to model spectral absorption, transmission, emission and scattering characteristics of the atmosphere (see reference [34] and references therein). The standard output of MODTRAN 4 consists of spectrally resolved transmittances, optical depths, radiances, and fluxes split into their components (e.g. transmittance of various constituents, direct and scattered radiation).

The paper is structured as follows: next section provides an overview of the radiative transfer governing equations and Monte Carlo solution methods. Section 3 describes the MOCRA code, including a presentation of the input quantities required to run a simulation and of the produced output, both within and outside the DOAS context. Sample calculations are presented in section 4 and validated through comparison with MODTRAN4 and observations. The summary of the results obtained and the perspectives for future applications of MOCRA, especially concerning the DOAS frame, are given in the conclusions section.

2. Monte Carlo solution of the radiative transfer equation

The radiative transfer equation (RTE) for time independent problems can be written in the integral form as in Eq. (1) [35]

\begin{matrix} I (r, Ω) = Γ (r_{s}, Ω) \exp [- \int_{0}^{| r - r_{s} |} d s^{″} k (r - s^{″} Ω)] + \\ \int_{0}^{| r - r_{s} |} d s^{'} Q (r - s^{'} Ω, Ω) \exp [- \int_{0}^{s^{'}} d s^{″} k (r - s^{″} Ω)] \end{matrix}

where

I (r, Ω)

is the specific intensity of radiation (or radiance) received in the direction Ω at the observer point r, defined as the product of the photons distribution function by the quantity

c h ν

_, with c the velocity of light in a vacuum,

h

the Planck’s constant and

ν

the radiation frequency.

Γ (r_{s}, Ω)

is the specific intensity of the radiation source,

k

is the extinction coefficient, accounting for both scattering and absorption contributions.

Q

is the source term defined by Eq. (2)

Q (r, Ω) = S (r) + \int_{4 π} d Ω^{'} k_{s} (r, Ω \cdot Ω^{'}) I (r, Ω^{'}) l .

The quantity $S$ in Eq. (2) can be defined as $S (r) = {k^{'}}_{a} B$ ; where $B$ is the Planck’s function [30], and $k_{a}^{'}$ is derived from the absorption coefficient $k_{a}$ by means of Eq. (3)

{k^{'}}_{a} = k_{a} [1 - \exp (- h ν / K T)]

where K is the Boltzmann constant and T is the absolute temperature.

A more detailed description of the radiative transfer equation can be found in references [36,37] which define the Stokes parameters to describe the effects of polarization.

As it can be seen in Eq. (1), the specific intensity of radiation is determined by two contributions. The first one is the specific intensity hitting the system’s surface at point $r_{s}$ . The second term accounts for emitted and scattered photons from each path element $d s^{'}$ between $r_{s}$ and $r$ along $Ω$ direction. Each term is exponentially attenuated by scattering and absorption along the path between the source point and the observer location. The photonic interpretation of radiation allows applying the models usually adopted for simulating particles transport to solve the RTE in the atmosphere; many such models rely on the Monte Carlo approach. The Monte Carlo modeling involves a free optical path (or optical depth) to be first simulated. The optical depth is the exponent of the attenuation factors in the RTE (Eq. (1)):

τ (r, r_{s}) = \int_{0}^{| r - r_{s} |} d s^{″} k (r - s^{″} Ω)

The Monte Carlo approach to the solution of atmospheric radiative transfer consists of simulating the trajectories and occurring events, that together we call “histories”, of each individual photon emitted by the radiation source. Being the histories independent of one another by definition, the Monte Carlo method allows for the straightforward parallelization and the optimization of the computing resources.

The adjoint transport equation is very similar to Eq. (1). The source terms $Γ (r_{s}, Ω)$ and $S (r)$ , in fact, are evaluated at the actual position of the source, but in the integration intervals the quantity $r_{s}$ is replaced by $r_{d}$ , the detector position, and the incoming and outgoing directions $Ω$ and $Ω^{'}$ are exchanged each other. Due to isotropy of the medium, however, the scattering kernel depends on these directions only through their scalar product, so that exchanging the order of directions does not change the value of scattering kernel.

Though from a macroscopic point of view the attenuation of a radiation beam can be described by an analytical equation (Eq. (1), a single process occurring to a single photon has to be treated in a probabilistic way and classified as a random process. In this sense, every process a photon undergoes in the atmosphere can be described in terms of probability distribution functions (pdfs) [38]. The Monte Carlo approach consists in simulating the physical system by random sampling from the pdfs describing each specific process. Random samples from arbitrary pdfs can be obtained in different ways [39], but particular attention should be paid in using appropriate random number generators to avoid long-range correlations, especially in parallel Monte Carlo codes [40,41]. One of most common methods is sampling via inversion of the cumulative distribution function (cdf). Since the cdf, F(x), of an arbitrary pdf, f(x), is always uniformly distributed in the interval [0, 1], it is possible (i) to sample a random number R uniformly distributed between [0, 1], (ii) to equate the random number with the cdf (F(x) = R), and (iii) to invert the cdf and solve for x (x = F⁻¹ (R)). This paper uses the aforementioned method to randomly sample the quantities of interest, such as the distance between two consecutive collisions, the scattering angles (i.e., the direction after collision), the absorption rates, etc., from the various pdfs that describe the photon interactions within the atmosphere. The formal treatment of this procedure, here omitted for brevity, is described in references [42,43].

3. Description of MOCRA

MOCRA is a FORTRAN code which, given a monochromatic electro-magnetic radiation in the UV-near Infrared frequency range, performs a Monte Carlo simulation of the radiative transfer according to the physical and geometrical data of the assigned atmospheric environment. It is suitable to perform simulations of passive remote sensing measurements, as for instance DOAS observations. In this case a backward Monte Carlo scheme is adopted to simulate the UV-Visible radiance collected by a detector. The code simulates the atmosphere as a spherical multi-layer environment, taking into account also horizontal variations of physical parameters through three-dimensional multi-region rectangular and spherical geometries. The surface albedo is taken into account, assuming a Lambertian surface. Refraction is simulated according to Snell’s law.

In order to improve the performance of MOCRA, semi-analytical estimates of specific quantities (e.g., the absorption rates) and variance reducing techniques such as forced collision and expected leakage, are foreseen.

To evaluate the effects of small perturbations in the physical components of the atmosphere, a reference environment and several perturbed scenarios can simultaneously be simulated in a single run of the code. For each environment of interest in a simulation, the external library of physical data (interaction coefficients, phase functions and refractivity indexes) is required.

In the following subsections the main features of the MOCRA code are described: section 3.1 describes the radiation source and detectors representation; in section 3.2 the geometric and physical treatment of the atmosphere is depicted; section 3.3 concerns the photon history tracking, the end of history options and the variance reducing techniques; section 3.4 deals with perturbative analysis and section 3.5 gives an outlook of the quantities calculated by the model.

3.1 Radiation source and detectors

MOCRA assumes the Sun as the monochromatic source of radiation. The incoming photon directions are given by means of direction cosines (u_s, v_s, w_s), which can be either assigned by the user, or calculated on the basis of the astronomical parameters [44]. Through a basic knowledge of the standard meridian longitude with respect to the Greenwich meridian together with the longitude and latitude of the observer, the solar zenith and azimuth angles can be computed once the day of the year and the standard solar time are assigned. In this regard, the observer local time can be either assigned by the user or computed on the basis of the standard solar time, standard longitudes and the equation of time, accounting for the non-uniform terrestrial motion.

The radiation receivers handled by MOCRA can be schematized as point or disk detectors. A special treatment is foreseen for upward looking point detectors at ground level, in order to better take into account the refraction phenomenon at twilight when performing single scattering radiance calculations. In the following, the main receiver types allowed by MOCRA and the way in which collection of solar radiation is formally treated are described.

Vertical upward looking point detectors. Let I(z) be the radiance intensity due to a solar ray reaching the detector from the z-height level. Integrating the contributions from all z levels from the ground to the top boundary of the atmosphere gives the total intensity I_T, as in Eq. (5):

I_{T} = \int I (z) d z

or equivalently in Eq. (6):

I_{T} = \sum_{i} Δ z_{i} \int_{Δ z_{i}} I (z) (d z / Δ z_{i})

where the entire vertical range is subdivided into N intervals of thickness

Δ z_{i}

(i = 1,…,N) and the integrand is changed accordingly.

In the Monte Carlo simulation, a stratified sampling procedure can be performed, that consists in randomly choosing a set of z points uniformly distributed within each height subinterval $Δ z_{i}$ , and in estimating the corresponding $I (z)$ value. The average intensity value over the points belonging to the i-th subinterval times the vertical layer thickness $Δ z_{i}$ gives the contribution of that layer to the total intensity, the latter accounting for both single and multiple scattering radiance terms. Single scattering radiance ( $I_{S}$ ) is calculated using a forward Monte Carlo simulation of the photon path from the solar source coordinate to the altitude of each selected arrival point on the z-axis taking into account refraction. The effect of refraction is particularly important for SZA near or greater than 90° (horizon) and for solar rays crossing the lower atmospheric layers. Multiple scattering radiance ( $I_{M}$ ) is calculated using a backward Monte Carlo simulation from the selected vertical points to the sun. The albedo phenomenon is foreseen assuming a Lambertian surface, refraction is not accounted for along the path from the last scattering point to the sun.

Point detectors with assigned position and line of sight. In a O(x,y,z) reference system of coordinates, let P_d = (x_d, y_d, z_d) be the point detector spatial coordinates and (u_d, v_d, w_d) the direction cosines of its line of sight. Starting from P_d, a backward Monte Carlo simulation of a photon history with initial flight direction (u_d, v_d, w_d) is performed, taking into account both single and total scattering radiance.

Disk detectors with assigned diameter, Field Of View (FOV), position and orientation. Let P_d = (x_d, y_d, z_d) be the centre of the disk and (u_d, v_d, w_d) the direction cosines of the outward direction normal to the disk surface. By assuming D (≥ 0) the diameter of the disk and α (>0) its FOV, a fictitious cone with semi-amplitude α and the disk as a basis is generated. The vertex point P_v, ≡ (x_v, y_v, z_v) (possibly located outside the geometrical system) is assumed to be the virtual source point in a backward Monte Carlo simulation, with initial motion direction Ω for each photon uniformly distributed inside the cone. For each Ω direction, the crossing point with the detector disk gives the “true” starting point for the photon. Let ΔΩ be the solid angle amplitude and I(Ω) the contribution to radiance intensity collected by the receiver from the Ω direction. The total intensity I_T over all possible directions is assumed to be:

\begin{matrix} I_{T} = \int_{Δ Ω} I (Ω) d Ω \\ = Δ Ω \int_{Δ Ω} I (Ω) d Ω / Δ Ω \end{matrix}

Equation (7) can be simulated through a standard backward Monte Carlo estimating procedure for I_T. If D = 0, analogous considerations can be made, with the true starting point coinciding with the virtual one. If α = 0, the starting point is chosen uniformly on the disk and the initial direction coincides with the detector axis.

3.2 Geometrical and physical description of atmospheric environments

The atmosphere is modeled as a set of either (chemically) homogeneous or non homogeneous horizontal layers. In both cases, the albedo phenomenon is accounted for whenever the ground surface is hit. In case of non homogeneous horizontal layers, refraction is foreseen whenever the boundaries between contiguous regions are crossed. The physical characteristics of each region are provided to MOCRA with a data set contained in an external library. The various kinds of geometrical environments that are foreseen by the code are described in the following.

1) Horizontally homogeneous atmosphere. Starting from the ground, the atmosphere is divided into a series of horizontally homogeneous spherical shells, each one with its own physical properties. Two reference systems are simultaneously considered for the photon history tracking: O(x,y,z) and (R,φ), where R is the radius from the collision point to the earth centre and φ is the angle between the direction of vector $\vec{R}$ and the motion direction. The system based on (R,φ) coordinates is more suitable with respect to the Cartesian system to compute the escape probability and the leakage optical path for the travelling photon.
2) Horizontally non homogeneous atmosphere. Horizontal variations of the atmospheric parameters and different topography scenarios are simulated by the code using a multi-region rectangular or spherical geometry:
- a. Multi-region rectangular geometry. Given a rectangular domain (x_min, x_max, y_min, y_max), a set of sub-regions is generated by dividing the interval (x_min, x_max) in a number of subintervals. For each of them a subdivision of the interval (y_min, y_max) is performed. Each such generated sub-region has its own atmospheric environment and surface characteristics.
- b. Multi-region spherical geometry (Fig. 1 ). Let O(x,y,z) be a reference system with origin on the Earth centre and the z-coordinate axis oriented according to the zenith of the observer on the Earth surface. This geometry model foresees a set of “N” cones, each one with the vertex located at the centre of the Earth, the axis coincident with the z-coordinate axis, and the bisector of the vertex angle θ (< 90°) coincident with the zenith angle with respect to the vertical axis (Fig. 1a). It is assumed that θ₁ < … < θ_N. A set of “M” vertical half planes can be associated to each cone, except that to the first one: the cones axis (the z-coordinate axis) is the intersection of all half planes. As shown in Fig. 1b, for each vertical half plane, starting from the x-axis in the (x,y) coordinate plane, the azimuth angle ψ_j (0 ≤ ψ₁ < … < ψ_M = ψ₁ + 360°) is defined, being ψ_{k + 1} - ψ_k < 180° (k = 1, …, M - 1). In this way, the i-th geometrical space (i > 1) between the i-th and (i + 1)-th contiguous cones is subdivided in a number of (M_i - 1) regions. Figure 1c shows the horizontal section of a possible geometry with 3 cones and 5, 7 half planes for the second and third cone respectively. Each region generated according to the previous procedure, has its own vertical subdivision as previously described for the spherical multi-layer geometry, so that different ground levels can belong to different zones. Figure 1d shows a vertical section of conic regions with bisectors of vertex angles θ₁, θ₂, θ₃ and ground heights h₁, h₂ and h₃. A simple application of this type of geometry in MOCRA has already been described in reference [45].
  
  Fig. 1 Multi-region 3D spherical geometry: a) a set of three cones with the vertex at the Erth centre; b) vertical half planes having the z axis as intersection line with azimuth angle ψi and ψj defined in (x,y) plane; c) a horizontal section of a possible geometry with three cones and 5, 7 half planes for the second and third cone, respectively; d) a vertical section with three cones of aperture θ1, θ2, θ3 and ground heights h1, h2 and h3.
  Download Full Size | PDF

Whether the atmosphere is schematized as a horizontally homogeneous spherical multi-layer or a horizontally non homogeneous multi-region pattern, its vertical structure is described by the vertical profile of up to fifteen molecules and up to four distinct aerosol types, depending both on different size distribution and number density and on material type related to the different types of environment, and with various kinds of clouds and/or fogs. This information is supplied by an external library of physical data that describes the vertical composition of the atmosphere and provides the fundamental physical constants and aerosol scattering phase functions. Also the aerosol type is defined by the library. In the photon diffusion process, an essential role is played by molecular Rayleigh scattering, which affects the distance travelled by the photon and its motion direction following a collision. A theoretical formula of the scattering coefficient $k_{M s}$ which takes into account the dependence on wavelength $λ$ , refraction index and molecular number density is given by Kondratyev [46] together with phase function and refraction index expressions. In the MOCRA code, the following formula is employed, the same used in LOWTRAN-MODTRAN series codes [47]:

k_{M s} (λ) = \frac{24 π^{3}}{N λ^{4}} {[\frac{n^{2} - 1}{n^{2} + 2}]}^{2} \frac{6 + 3 δ}{6 - 7 δ}

where

δ

is the depolarization factor, assumed equal to 0.0279 as in the MODTRAN codes. It is apparent that Rayleigh scattering coefficient rapidly increases as the radiation wavelength decreases. The following formula is used for the phase function [48]

p (θ) = \frac{3}{16 π} \frac{2}{(2 + δ)} [(1 + δ) + (1 - δ) \cos^{2} θ]

which is directly computed by MOCRA [46]. With reference to the aerosol reaction coefficients and phase functions, the Mie theory is applied, which holds for spherical homogeneous particles with radii

r

comparable to the wavelength

λ

of the incident radiation, but, being this theory very complex from both analytical and numerical points of view, option is left to use the well known Henyey-Greenstein distribution

p (θ) = \frac{1}{4 π} \frac{1 - g^{2}}{{(1 + g^{2} + 2 g \cos θ)}^{\frac{3}{2}}}

on the basis of the asymmetry factor g, given by

g = \iint_{4 π} p (θ) \cos θ d Ω

the value of which is given in the library. As an option, it is possible to exclude the effects of either molecular or aerosol particles in a simulation.

3.3 Path tracking and variance reducing techniques

History tracking in Monte Carlo simulations is based on sampling random variables distributed according to appropriate, discrete or continuous, pdfs which describe the random walk of photons and their interactions with atmospheric constituents. As mentioned above, the present work uses sampling via inversion of cdfs.

To trace, in its fundamental lines, a photon history, the initial weight (w₀), position and motion direction have to be assigned to the starting photon. The statistical weight carries information about both intervening events on the photon itself (in general the photon weight is reduced at each interaction) and about the computational techniques or variance-reducing strategies adopted in the simulation. At the beginning of the history and after each collision the photon undergoes with an atmospheric constituent, the optical path to be traveled before the next interaction is calculated, and the geometrical path length is evaluated accordingly from the knowledge of the interaction coefficient of the crossed layers. Using a geometry routine, it is decided whether the photon exits a layer or crosses a vertical boundary separating horizontally different regions or it undergoes an interaction within the layer. In case of leakage from the atmosphere (i.e. the photon crosses the last layer or the Top of the Atmosphere), the current photon weight is stored in a suitable counter and end of history calculations are performed. When an interaction occurs at a given point, a new weight for the photon is computed as follows: $w_{n} = w_{n - 1} \frac{(k_{A s} + k_{M s})}{k}$ with $k_{A s}$ , $k_{M s}$ and $k$ the scattering coefficients for aerosols and molecules and the total interaction coefficient, respectively, at the actual collision point, being $w_{n - 1}$ the photon weight before collision. A weight value proportional to the absorption probability is stored in a suitable estimator, thus obtaining a semi-analytical estimate of absorption. This is the so-called “collision estimator” for absorption which is updated every time a collision occurs. An estimate of the absorption can also be performed by means of a “distance estimator”, that is, through the computation of the product between the photon path length in a given layer and the mean absorption coefficient along that path. Since this estimation makes use of the average absorption cross section of the considered layer rather than its punctual value, this estimator gives an approximated value of the absorption contribution. On the other hand, it has the advantage of supplying analytical estimate of absorption independently of collision to occur. According to the scattering probability due to an aerosol particle $p_{A s} = \frac{k_{A s}}{k_{A s} + k_{M s}}$ and to a molecular particle, $p_{M s} = \frac{k_{M s}}{k_{A s} + k_{M s}}$ , the choice between aerosol and molecular scattering is performed. In case of aerosol scattering, the kind of aerosol involved is selected among 4 types (allowing to include, together with the main rural, urban or maritime aerosol, other aerosol coming from specific sources such as smokestacks or fires) based on their probabilities $p_{A i s} = \frac{k_{A i s}}{k_{A s}}$ , with i = 1 to 4. On the basis of molecular or aerosol scattering phase functions, a new motion direction is selected and a new path starts. If the distance $d$ to be travelled is greater than the distance from the new starting point to the boundary of the actual layer along the motion direction and a leakage phenomenon from the atmosphere does not occur, the optical path exhausted in the layer will be computed, the remaining one being utilized for the next path in the new layer. If refraction is foreseen, according to the refraction indices $n_{1}$ and $n_{2}$ belonging to the two contiguous layers, the new refracted motion direction is computed except that when $n_{1} > n_{2}$ and the angle of the incident direction on the boundary between the two layers is greater than $θ_{L} = arc \sin \frac{n_{2}}{n_{1}}$ . In this case a total reflection occurs and the path goes back into the previous layer with direction defined by the Descartes’ law. If the photon hits the ground surface and the albedo phenomenon is foreseen in the calculation with an albedo coefficient α, a photon with weight $w = α w$ is reflected and the new motion direction is chosen randomly from the Lambertian distribution function.

The history tracking continues until one among possible user selectable end of history criteria, each one depending on the type of problem to be solved and the degree of precision required, is satisfied. Three end of history options are foreseen in MOCRA. 1) Given a minimum weight fraction, $w_{\min}$ , the history will end if $\frac{w}{w_{0}} < w_{\min}$ . The minimum weight fraction has to be carefully chosen to avoid expensive calculations giving inessential contributions to the estimator. A value of $w_{\min}$ equal to 10⁻⁵ is generally adopted. The sum over the histories of the unprocessed weights will give a measure of the bias introduced adopting an arbitrarily chosen minimum weight. 2) Given a cutting weight value, $w_{c u t}$ , if $\frac{w}{w_{0}} < w_{c u t}$ after a collision, the photon weight will no more be changed. When a next interaction occurs, if it is an absorption event the history tracking will end, otherwise it will continue without further weight reduction. A $w_{c u t}$ value equal to 0.2 is generally assigned. 3) Given a cutting weight value $w_{c u t}$ , if it is $\frac{w}{w_{0}} < w_{c u t}$ a “russian roulette” game is played: with a probability $p_{r} = \frac{w}{w_{0}}$ the photon will survive and the initial weight $w = w_{0}$ will be restored; with probability $1 - p_{r}$ the photon will be killed and the photon weight $w$ will be lost. At the end of all processed histories, the sum of the total weight gained should be statistically equal to the total weight lost. In this case $w_{c u t} = 0.5$ is a reasonable value that avoids high jumps in the weight, which can affect the variance.

Batches of pre-assigned number of photon histories are generally run in a MOCRA simulation; each batch is assumed as a single numerical experiment. Besides being more appropriate than a simulation of all photon histories in a unique model run as regard the variance calculation, this procedure allows high number of histories (as it should be necessary) to be run without significant precision loss when collecting the results coming from the simulation [49]. A restart option is also foreseen, to make the results of an assigned number of batches available for a subsequent continuation of the simulation. To make the comparison between slightly different calculations (in the physical or geometrical parameters) sufficiently robust, a strategy on the random number generation is adopted for both batches and histories inside a batch so that, in each calculation, corresponding histories begin with the same random number. This guarantees that unaltered corresponding histories in the two calculations give the same results for the same quantity.

Analytical estimates of specific quantities can be performed in the simulation to reduce as much as possible the variance of the calculation. Absorption rates, for instance, are estimated semi-analytically through the distance estimator described above. Another example of analytical estimators is the leakage estimator which is evaluated when a collision occurs. At each collision point the optical distance $τ_{e}$ from that point to the outer atmospheric boundary along the motion direction is computed. The escaping weight $w_{e} = w e^{- τ_{e}}$ is the analytical contribution to leakage from the collision point. The sum over all collisions of these contributions gives the analytical estimate of leakage for that history. The colliding weight $w_{c} = w (1 - e^{- τ_{e}})$ is assigned to the photon and used for evaluating the collision estimator in the following part of the history.

It should be noticed that, despite the variance is a function of the number of histories, N, its reduction cannot be achieved simply increasing N, because this would increase the computational time. To reduce the variance of a factor 1/F, in fact, it is necessary to multiply the number of histories by a factor F² [38]. Therefore, in addition to analytical evaluations, that certainly reduce the computing time, variance reducing techniques (vrt) are adopted. They have the property of somewhat altering the natural sequence of the events a photon would undergo, in order to increase the photon statistics saving also computing time, giving rise to a more efficient contribution to the required estimates. Probably the most common vrt, also used in MOCRA, is forced collision: for motion directions that do not hit a reflecting boundary, the leakage optical distance $τ_{e}$ is computed and the probabilities for a photon to escape from the system ( $e^{- τ_{e}}$ ) and to collide inside it ( $1 - e^{- τ_{e}}$ ) are computed. In this way, a travelling photon with weight $w_{e} = w e^{- τ_{e}}$ will escape from the system, while a travelling photon with weight $w_{c} = w (1 - e^{- τ_{e}})$ will collide. It is worth pointing out that both photons fates are considered when the absorption distance estimator is calculated. When a forced collision is simulated, the optical path before collision is randomly selected from the pdf $f (x) = \frac{e^{- x}}{1 - e^{- τ_{e}}}$ (truncated exponential distribution) rather than from the simply exponential distribution, and the photon travels with weight $w = w_{c}$ . End of history criterion 1 is adopted in this case and, as a precaution, a maximum number of collisions per history is either assigned by the user or assumed equal to 1000. Forced collision, however, is a time consuming technique and may cause considerable changes in the travelling photon weight which can negatively affect the variance of the calculation. Therefore, it should only be used as far as information for a specific layer (or region) of the atmosphere is not correctly achievable using a standard simulation (for layers with small optical thickness, for instance, for which an unbiased simulation would give rise to small interaction probabilities). A Figure Of Merit (FOM) can be defined to simply measure the efficiency of the simulation process:

F O M = \frac{1}{t σ_{μ}^{2}}

where

t

is the computing time and

σ_{μ}

is the standard deviation of the estimated parameter.

A more comprehensive review of variance reducing techniques can be found in references [28,50].

3.4 Perturbative calculations

Small perturbations in the physical atmospheric components, with respect to a reference atmosphere, and their effect on radiation transport, can be accounted for in MOCRA. To do so, photon histories are traced sampling the intervening events from the pdfs describing both the unperturbed and the perturbed situation. To this aim an “importance sampling” Monte Carlo technique [42], which allows defining proper weight factors for the photons travelling in the perturbed and unperturbed environments is employed. In the radiation transport simulation, let $k$ and $k^{'}$ be the extinction coefficients in a given region for the unperturbed and the perturbed environment, respectively. When crossing a boundary of an atmospheric region in which the perturbative calculation is performed, a weight factor equal to $\exp (- (k^{'} - k) d)$ has to be associated to travelling particle in the perturbed environment, being $d$ the distance from the starting point to the boundary, since the probability to leave the region is given by $\exp (- k^{'} d)$ in the perturbed environment and $\exp (- k d)$ in the unperturbed one. If a collision occurs in the region, the weight factor will be $(k^{'} / k) \exp (- (k^{'} - k) d)$ , where $d$ is the distance traveled in the region. $k^{'} \exp (- k^{'} d)$ is the probability for the photon to collide after a distance $d$ is travelled in the perturbed environment and $k \exp (- k d)$ is the analogous quantity in the unperturbed environment. Therefore, when perturbative calculations are carried out, at least 2 weights are attached to each photon during a simulation, one for each situation involved.

Estimators of the differential effects of the perturbation on certain quantities are computed by collecting, during the history tracking, the differences in the quantities of interest between the perturbed and unperturbed environments, so that the effects of small perturbations, difficult to evaluate with separate calculations, are computed in a single run of the code. The maximum number of perturbed environments that can be handled by the code is limited only by memory usage. It must be remarked that, according to this technique, the leading photon history is governed by the unperturbed system, so that an end of history in this system will cause the end of histories also in the perturbed environments. Perturbations are not foreseen to affect the angular distributions. When horizontally non homogeneous atmospheric conditions are considered and perturbations are foreseen in different regions, the unperturbed environment must be the same for all of them.

3.5 Computed quantities

The first quantities computed by the MOCRA code are the atmospheric optical depth $τ$ , defined by Eq. (4), the atmospheric transmittance, defined as $T r = e^{- τ}$ , and the spatial distribution of single and total scattering irradiance reaching the ground surface. Besides that, the single and total scattering radiance and the ground reflected radiance, received by a simulated detector are evaluated. Option is left to avoid scattered radiance or reflected radiance calculation. The aerosol and molecular contributions to scattering are separately accounted for. The spatial distribution of the single and total scattering contribution to the radiance is calculated along the detector axis, according to the geometrical subdivision assigned by the user along the detector axis. Radiance calculations using a source other than the sun, placed at an assigned distance from a point detector, can be performed. For each molecules of interest, absorption vertical optical depths, intensity weighted optical paths for the single and total scattering radiance and for the reflected radiance, and the AMF are calculated. The AMF of a specific trace gas, in particular, is used to convert the SCD measured with DOAS instruments into a VCD. Several types of RTMs have been developed for such an aim taking into account either single (e.g., the RTM AMEFCO [23]), or multiple (e.g., the PROMSAR RTM [26]) scattering. Being I and I* the radiance (single or total scattering) detected by a receiver with and without the trace gas of interest in the atmosphere, respectively, the AMF of the species can be calculated as in Eq. (13) [51]:

A M F = - \frac{\ln (\frac{I}{I^{*}})}{δ_{a}}

where δ_a is the vertical absorption optical depth of the molecular species of interest. Equation (13) is a general expression for the AMF of a given trace gas, which holds especially in case of strong absorbers and will be thus referred to also as AMF_SA. For weak absorbers an approximating formula for the AMF (AMF_WA) is given by Eq. (14) (see also reference [51]):

A M F_{W A} = \frac{δ_{O P T}}{δ_{a}}

where δ_OPT is the intensity weighted absorption optical path from the sun to the detector defined as

δ_{O P T} = \frac{\sum_{i} D_{i} I_{i}^{*}}{\sum_{i} I_{i}^{*}}

where

D_{i}

is the distance traveled by photons in the i-th layer,

I_{i}^{*}

is the radiance collected by the detector from the i-th layer in the absence of the examined gas.

Due to the possibility for the MOCRA code to perform simultaneously radiative calculations with and without the molecular particle of interest (perturbative calculation), the AMF can be calculated as in Eq. (13). The approximated value of the AMF can altogether be calculated (Eq. (14), and compared with the former to check the hypothesis of weak absorption. An example of such a comparison is shown in Fig. 2 for the O₃ and NO₂ AMFs at 310 nm, for a ground based vertically looking detector.

Fig. 2 AMFs for O3 (left) and for NO2 (right) calculated as in Eq. (13) (AMFSA, red line) and Eq. (14) (AMFWA, blue line). The red line for NO2 has been shifted upward of 0.5 to make it distinguishable from the blue line.

Download Full Size | PDF

As expected, significant differences in the AMF of O3 (left panel) between the two formulas can be observed, due to the fact that ozone is a strong absorber in the near UV. The difference between the exact and approximated AMF is even greater than 50% for solar zenith angles (SZA) near or greater than 90°. On the other hand, NO₂ is a much weaker absorber than ozone at 310 nm, leading to differences always less than 1% (0.1%) between the two AMF calculations in case of single (multiple) scattering simulations (the red line for NO2 AMF in the figure is shifted of 0.5 upward).

A very useful, quite recent concept is that of box-AMF [27]. The box-AMF for a certain atmospheric “box” can be calculated as the derivative of the SCD within that box with respect to the VCD [21]. This is equivalent to calculate the derivative of the logarithm of the intensity with respect to the number density n_b of the gas in the box or, more precisely, the absorption coefficient β_b, as indicated in Eq. (16).

A M F_{b} = \frac{d S C D_{g}}{d V C D_{b}} = - \frac{1}{h_{b} σ_{b}} \frac{d \ln I_{g}}{d n_{b}} = \frac{1}{h_{b}} \frac{d \ln I_{g}}{d β_{b}} .

It can be noted that Eq. (16) reduces to Eq. (13), provided that $δ_{a}$ is substituted by its value $δ_{a} = h_{b} β_{b}$ for the box. As for the AMF, the MOCRA code calculates the box-AMF and its approximated expression for weak absorbers (or little-explored regions), which gives a reasonable estimate of this factor with less computational demand. The approximated box-AMF for a specific trace gas is given by Eq. (14), provided that $δ_{a}$ and $δ_{O P T}$ are calculated for the box, so that $δ_{a} = h_{b} β_{b}$ and in Eq. (15) $D_{i}$ is supposed equal to zero if the layer is outside the box.

When perturbative calculations are not performed, the approximated expression for the box-AMF can still be applied, provided that $I_{i}^{*}$ in Eq. (15) is replaced by $I_{i}$ , i.e. the radiance detected in the presence of the gas of interest. Several boxes can be considered in a single simulation. The assumption of monochromatic photons and the energy dependence of absorption are probably the reasons why the ozone AMF in the visible wavelength range (the so-called “Chappuis band”, see section 4) presents a strong variability with wavelength. This gives rise to some problems in using these results for the interpretation of passive DOAS measurements, which would be the subject of future investigations.

4. Results, applications and comparison exercises

This section gives an outlook of possible applications of the MOCRA model, especially in the frame of DOAS. However, to verify the reliability of the code in solving radiance problems, two comparison exercises between MOCRA and MODTRAN4 for a variety of reciprocal positions between the instrument line of sight (u_d, v_d, w_d) and the solar direction (u_s, v_s, w_s) and different atmospheric visibility values, are first presented. Some applications of MOCRA in the frame of DOAS will be subsequently shown, for which a deterministic model, like MODTRAN4, should not obtain sufficiently accurate results due to multiple scattering effects and geometric requirements. The first exercise consists of simulating a point detector, located at an altitude of 1 km above the ground (as if it was flying on board a research aircraft) with lines of sight $α = 180 °$ and $α = 150 °$ (with respect to the vertical upward direction). The solar zenith angles are $θ = 30 °$ and $θ = 60 °$ , while azimuth angles are $ψ = 0 °$ and $ψ = 90 °$ between the vertical plane containing the detector line of sight and the sun direction. In Tables 1 and 2 the obtained results are quoted for a visibility of 5 and 23 km, respectively. The considered atmospheric scenario chosen for this simulation is that of MODTRAN4 for mid-latitude winter conditions, for a wavelength of 550 nm. The solar constant for that wavelength is calculated with MODTRAN4 as 6.1134 10⁻⁶ Watt/(cm² cm⁻¹), a typical value for January. The mean surface albedo (0.3) was assumed.

Table 1. Point Detector: MOCRA - MODTRAN-4 Comparison for vis = 5 Km (for wavelength and atmospheric conditions see text)*

View Table | View all tables in this article

Table 2. Point Detector: MOCRA - MODTRAN-4 Comparison for vis = 23 Km^a

View Table | View all tables in this article

For each combination of $α, ψ, θ$ and the two visibility values, the compared radiometric quantities are single scattering and total scattering radiance, reflected radiance at ground, along with the sum of total and reflected radiances. For each quantity the percent difference between the MOCRA and MOTRAN4 estimates is reported. Discussing the details of the comparison as a function of the different values of the input parameters is outside the scope of the paper. It is important to stress that a satisfactory agreement between MOCRA and MODTRAN4 for both scattered and reflected radiances has been found. The difference between the two codes is less than 1% for single scattering and less than 3% for multiple scattering. Monte Carlo methods have proven to be particularly suitable to simulate radiative transfer problems with marked multiple scattering anisotropy [42] and geometrically confined sources with narrow fields of view, that can give rise to troubles known as “ray effect” and “false scattering” [52], if not properly treated.

A second comparison exercise between MOCRA and MODTRAN4 was performed assuming a horizontal looking ground-based point detector. In this comparison exercise we focus on single scattering calculations only. In this way, the effects of using different methods to solve the RTE when multiple scattering is included are not accounted for and do not contribute to possible discrepancies in the results. A mid-latitude winter model for the atmosphere was used in this exercise. The considered radiation wavelengths were 310 nm, 435 nm, 550 nm and 700 nm. Seven visibility values, chosen by considerations about the relationship between visibility and aerosol load, (5 km, 7.6 km, 9.5 km, 13 km, 18 km, 23 km and 35 km), solar zenith angles of 30° and 60°, and solar azimuth angles of 0° and 90° have been considered. Being too long to show all results, they are summarized in Table 3 as the average values over all visibilities, solar zenith angles and azimuth angles. The table shows, for each radiation wavelength and the respective value of the solar constant calculated with MODTRAN4, the average radiance values computed with MOCRA and MODTRAN4, along with their mean and standard deviation differences. The mean difference between the two codes is relatively low, less than 3%, with maximum standard deviation of about 5% at 310 nm, suggesting a really good agreement between the two models in calculating the single scattering radiance reaching such a simulated detector in a cloud-free atmosphere [25].

Table 3. Horizontal Looking Point Detector Radiance, MOCRA-MODTRAN 4 Comparison, Mid Latitude Winter Model, Urban Extinction. Units: watt/ (cm²-ster-cm⁻¹)

View Table | View all tables in this article

In what follows, we move to some practical uses of the MOCRA model with particular emphasis on DOAS applications: calculation of the AMF and simulation of measured SCDs.

Figure 3 shows the O₃ (left) and NO₂ (right) AMFs as a function of the SZA at 310 nm (blue), 435 nm (green), and 550 nm (red), calculated with MOCRA for a vertical looking detector at ground.

Fig. 3 AMFs computed for ozone and nitrogen dioxide at three different wavelengths: 310 nm (blue), 435 nm (green), 550 nm (red). The atmospheric model used was that of MODTRAN4 for mid latitude winter, with urban aerosol extinction, visibility = 35 Km.

Download Full Size | PDF

In the zenith-looking geometry, the AMF of stratospheric absorbers like ozone and NO2 increases with solar zenith angle for ground-based observations, due to increasing light path in the upper atmosphere. The AMF has a maximum (higher at 435 nm and 550nm than at 310 nm) near 90 ° SZA. The dependence of the AMF on wavelength is related to the fact that Rayleigh (e.g., molecular) scattering is a strong function of wavelength and also the absorption varies with wavelength. Figure 3 shows that at higher SZA (SZA greater than about 70°) there is a remarkable difference between the AMF of both gases at 310 nm (UV) and at 435-550 nm (visible), and the difference increases with increasing SZA. At low sun, in fact, the AMF is smaller in the UV than in the visible as more UV light is scattered before travelling the long distance in the atmosphere.

The choice between spectral intervals in which DOAS measurements of a specific gas are performed is usually driven by considerations about the absorption cross section of the gas. It is obviously important to perform measurements of a gas in a spectral range in which its absorption cross section is quite high (strong absorption). This requirement is not sufficient, however: it is also required that the absorption features of the trace gases are easily recognizable by the dedicated elaboration software. Figure 4 shows the absolute absorption cross sections (ACS) of O₃ and NO₂, plotted together in the 250-650 nm spectral range. The absorption cross section of O₃ is higher in the near UV (Hartley-Huggins band) than in the visible (Chappuis band). Chappuis band, however, is usually chosen for ozone retrieval mainly because the O3 ACS is less temperature dependent and the AMF is less wavelength dependent here than in the UV. The spectral range around 430 nm is normally chosen for DOAS NO₂ measurements.

Fig. 4 Absolute absorption cross sections for O3 and NO2 in the 250-650 nm spectral range at a 273 °K temperature.

Download Full Size | PDF

The second example of application of the MOCRA model discussed in this paper is the simulation of the slant column densities of a given trace gas, and the comparison with real observations. Such a comparison is useful to make assumptions on the vertical profile of the trace gas whose SCD is simulated by the model, as well as the atmospheric conditions at the time of measurements (such as the aerosol load, the atmospheric visibility, the presence of clouds, etc., see e.g [25,51].). In the example reported on here, the modeled SCD for NO2 are compared to experimental values obtained at the observatory of the Geophysics centre of Evora (38°33′36”N, 7°54′ W), Portugal, on 22 August 2009 in the p.m. part of the day. The atmospheric scenario for the MOCRA simulation was that of mid-latitude summer conditions with the visibility value set according to the measurements performed by a nephelometer placed in the measurement site. Figure 5 shows the comparison of simulated (blue symbols) and measured (cyan symbols) NO₂ SCDs. The percent error is on average less than 20% for solar zenith angles less the 56°, and increases for greater solar zenith angles, though keeping the agreement between the model and observations reasonably good.

Fig. 5 Comparison between simulated (blue) and measured (cyan) NO₂ slant column densities. The measurements were carried out on 22 august 2009 at Evora (38°33′36”N, 7°54′ W), Portugal.

Download Full Size | PDF

5. Summary and conclusions

In this paper, the radiative transfer model MOCRA has been described. It is a FORTRAN code which performs a Monte Carlo simulation of radiative transfer in a very general three-dimensional atmospheric scenario. When simulating passive remote sensing measurements, such those performed by a DOAS spectrometer, the MOCRA model adopts a backward Monte Carlo scheme, which is a useful tool to overcome problems arising from the narrow field of view of the detector. The main advantages of this code can be resumed in the following items:

1) the possibility to take into account not only vertical, but also horizontal variations of the physical and chemical atmospheric parameters, along with topographic features of the surface, through the implementation of three dimensional multi region rectangular and spherical geometries;
2) the possibility to evaluate the effects of small atmospheric perturbations in a single model run (perturbative calculations). In this way, small differences in some key parameters that influence the radiation transport and photons interactions with the atmosphere can be accounted for. This has direct consequence when the model is used to interpret DOAS remote sensing observations, since the AMF (and box-AMF) of a target trace gas can be evaluated through the exact definition formula. A number of different atmospheric perturbations can be accounted for simultaneously in a single run simply defining, for each of them, appropriate perturbed environments with respect to a reference scenario.

Other interesting features of MOCRA are the realistic treatment of multiple scattering, which relies on the employed Monte Carlo approach, and the use of variance reducing techniques to obtain statistically significant results with reasonable computing times. The choice of the particular variance reducing technique to be implemented depends on the problem at hand and strongly relies on the user experience. Forced collision is mostly used in MOCRA.

In the present paper, MOCRA has been first compared to the state-of-the-art model MODTRAN4, giving overall a very satisfactory agreement in the radiative quantities we have focused on, such as the single and total scattering radiance and the ground reflected radiance, for a variety of atmospheric scenarios and measurement configurations. For the single and multiple scattering radiance, for instance, the percent differences between the two models are less than 1% and less than 3% respectively. Multiple scattering is differently implemented in MOCRA with respect to MODTRAN4, leading to higher, though still small, differences between the two models in the multiple than single scattering radiance calculations.

Results for DOAS-related quantities, the AMF and SCD of selected trace gases, obtained with MOCRA have also been discussed in the paper. The AMF is the primary quantity calculated with the RTMs to interpret DOAS measurements of diffuse solar radiation, while simulation of SCDs and comparison with actual measurements is useful to make hypotheses on the vertical profile of the target species at the time of measurements. In this paper we first discuss the difference between the exact and approximated formula of the AMF (or box-AMF) and check the hypothesis of, respectively, strong and weak absorption under which the two formulas are valid. Then we show that MOCRA is able to capture the dependence of the AMF on the wavelength of radiation for the specific case of stratospheric absorbers, like ozone and nitrogen dioxide, and zenith looking geometry. As an example of application of the model, we compare the SCDs of NO2 simulated with MOCRA with SCD observations performed at the observatory of the Geophysics Centre of Evora in the south of Portugal during a typical day of the 2009 summer season. The agreement found between the model and observations can be considered as satisfactory, with percent error of about 20% for SZA lower than 56° and increasing for greater SZA, taking account of the difficulty to include in the simulation atmospheric inputs as close as possible to the reality. The only input parameter for the model which was constrained by observations is the visibility, due to the availability of concurrent nephelometer measurements performed in the same site.

In the frame of DOAS, the MOCRA model provides the main input values required to retrieve trace gases VCDs, AMFs, box-AMFs and vertical profiles by means of appropriate algorithms applied to the SCDs directly retrieved from the measurements of diffuse solar radiation. The paper has discussed some examples of application of this model to typical DOAS configurations, simulating measurements that could be performed with ground-based zenith or horizontal looking detectors, or detectors placed at a certain atmospheric height with respect to the sea level, or directly comparing simulations with SCD observations.

Due to the 3D spherical geometry, moreover, MOCRA has already given preliminary interesting results in computing the NO2 box-AMFs for limb viewing measurements performed with the SCIAMACHY (Scanning Imaging Absorption Spectrometer for Atmospheric CHartographY) instrument onboard ENVISAT [53]. These preliminary results deserve further investigations for their important application to limb-viewing satellite measurements in target atmospheric regions like the (still under-sampled) Upper Troposphere-Lower Stratosphere (UT-LS), a challenge region for most satellite measurements.

Acknowledgments

This work was partially funded by ENVIREN (Environmental Regional Network) laboratory and TRIL-ICTP (Training in Italian Laboratories-International Centre of Theoretical Physics) Trieste.

The paper was also financially supported through FEDER (Programa Operacional Factores de Competitividade – COMPETE) and National funding through FCT – Fundação para a Ciência e a Tecnologia in the framework of project FCOMP-01-0124-FEDER-014024 (Refª. FCT PTDC/AAC-CLI/114031/2009).

The authors would like to thank E. Cupini for help in developing previous and present versions of the model and M. Ranieri for help in graphical realization of some figures of the paper.

References and links

1. J. Lenoble, ed., “Standard Procedures to Compute Atmospheric Radiative Transfer in a Scattering Atmosphere - Volume I, International Association of Meteorology and Atmospheric Physics (IAMAP), ” Boulder, Colorado, USA (1977).

2. Y. Fouquart, W. M. Irvine, and J. Lenoble, eds., “Standard Procedures to Compute Atmospheric Radiative Transfer in a Scattering Atmosphere - Volume II, International Association of Meteorology and Atmospheric Physics (IAMAP),” Boulder, Colorado, USA (1980).

3. K. Stamnes, S.-C. Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media,” Appl. Opt. 27(12), 2502–2509 (1988). [CrossRef] [PubMed]

4. G. N. Plass and G. W. Kattawar, “Monte carlo calculations of light scattering from clouds,” Appl. Opt. 7(3), 415–419 (1968). [CrossRef] [PubMed]

5. G. N. Plass and G. W. Kattawar, “Radiative transfer in an atmosphere-ocean system,” Appl. Opt. 8(2), 455–466 (1969). [CrossRef] [PubMed]

6. E. Tinet, S. Avrillier, and J. M. Tualle, “Fast semianalytical Monte Carlo simulation for time-resolved light propagation in turbid media,” J. Opt. Soc. Am. A 13(9), 1903–1915 (1996). [CrossRef]

7. A. De Matteis and R. Simonini, “A New Monte Carlo Approach to the Adjoint Boltzmann Equation,” Nucl. Sci. Eng. 65, 93–105 (1978).

8. A. De Matteis and R. Simonini, “A Monte Carlo Biasing Scheme for the Adjoint Photon Transport,” Nucl. Sci. Eng. 67, 309–316 (1978).

9. D. G. Collins, W. G. Blättner, M. B. Wells, and H. G. Horak, “Backward monte carlo calculations of the polarization characteristics of the radiation emerging from spherical-shell atmospheres,” Appl. Opt. 11(11), 2684–2696 (1972). [CrossRef] [PubMed]

10. W. G. Blättner, H. G. Horak, D. G. Collins, and M. B. Wells, “Monte carlo studies of the sky radiation at twilight,” Appl. Opt. 13(3), 534–547 (1974). [CrossRef] [PubMed]

11. B. Mayer, S. W. Hoch, and C. D. Whiteman, “Validating the MYSTIC three-dimensional radiative transfer model with observations from the complex topography of Arizona's Meteor Crater,” Atmos. Chem. Phys. Discuss. 10(5), 13373–13405 (2010), doi:. [CrossRef]

12. J. F. Noxon, “Nitrogen dioxide in the stratosphere and troposphere measured by ground-based absorption spectroscopy,” Science 189(4202), 547–549 (1975). [CrossRef] [PubMed]

13. U. Platt, D. Perner, and H. W. Patz, “Simultaneous measurement of atmospheric CH2O, O3 e NO2 by differential optical absorption,” J. Geophys. Res. 84(C10), 6329–6335 (1979). [CrossRef]

14. U. Platt and D. Perner, “Direct measurement of atmospheric CH2O, HNO2, O3, NO2, and SO2 by differential optical absorption in the near UV,” J. Geophys. Res. 85(C12), 7453–7458 (1980). [CrossRef]

15. H. K. Roscoe, P. V. Johnston, M. Van Roozendael, A. Richter, A. Sarkissian, J. Roscoe, K. E. Preston, J.-C. Lambert, C. Hermans, W. DeCuyper, S. Dzienus, T. Winterrath, J. Burrows, F. Goutail, J.-P. Pommereau, E. D’Almeida, J. Hottier, C. Coureul, R. Didier, I. Pundt, L. M. Bartlett, C. T. McElroy, J. E. Kerr, A. Elokhov, G. Giovanelli, F. Ravegnani, M. Premuda, I. Kostadinov, F. Erle, T. Wagner, K. Pfeilsticker, M. Kenntner, L. C. Marquard, M. Gil, O. Puentedura, M. Yela, D. W. Arlander, B. A. Kastad Hoiskar, C. W. Tellefsen, K. Karlsen Tornkvist, B. Heese, R. L. Jones, S. R. Aliwell, and R. A. Freshwater, “Slant column measurements of O3 and NO2 during NDSC intercomparison of zenith-sky UV-visible spectrometer in June 1996,” J. Atmos. Chem. 32(2), 281–314 (1999).

16. A. Petritoli, G. Giovanelli, I. Kostadinov, F. Ravegnani, D. Bortoli, P. Bonasoni, F. Evangelisti, U. Bonafe, and F. Calzolari, “Tropospheric and stratospheric NO2 amount deduced by slant column measurements at Mt. Cimone Station,” Adv. Space Res. 29(11), 1691–1695 (2002). [CrossRef]

17. L. M. Perliski and S. Solomon, “On the Evaluation of Air Mass Factors for Atmospheric Near-Ultraviolet and Visible Absorptiion Spectroscopy,” J. Geophys. Res. 98(D6), 10363–10374 (1993). [CrossRef]

18. S. Solomon, A. Schmeltekopf, and R. W. Sanders, “On the interpretation of zenith sky absorption measurements,” J. Geophys. Res. 92(D7), 8311–8319 (1987). [CrossRef]

19. J. Pukıte, S. Kühl, T. Deutschmann, W. Wilms-Grabe, C. Friedeburg, U. Platt, and T. Wagner, “Retrieval of stratospheric trace gases from SCIAMACHY limb measurements,” Proceedings of the First Atmospheric Science Conference, 8–12 May, ESA/ESRIN, Frascati, Italy, ESA SP-628, (2006). http://earth.esa.int/workshops/atmos2006/participants/1148/paper proc Frasc 2.pdf

20. S. Kühl, J. Puķīte, T. Deutschmann, U. Platt, and T. Wagner, “SCIAMACHY Limb Measurements of NO2, BrO and OClO, Retrieval of vertical profiles: Algorithm, first results, sensitivity and comparison studies,” Adv. Space Res. 42(10), 1747–1764 (2008). [CrossRef]

21. J. Puķīte, S. Kühl, T. Deutschmann, U. Platt, and T. Wagner, “Accounting for the effect of horizontal gradients in limb measurements of scattered sunlight,” Atmos. Chem. Phys. 8(12), 3045–3060 (2008). [CrossRef]

22. G. Hönninger, C. von Friedeburg, and U. Platt, “Multi axis differential optical absorption spectroscopy (MAX-DOAS),” Atmos. Chem. Phys. 4(1), 231–254 (2004), http://www.atmos-chem-phys.net/4/231/2004/acp-4-231-2004.pdf. [CrossRef]

23. A. Petritoli, F. Ravegnani, G. Giovanelli, D. Bortoli, U. Bonafè, I. Kostadinov, and A. Oulanovsky, “Off-axis measurements of atmospheric trace gases by use of an airborne ultraviolet-visible spectrometer,” Appl. Opt. 41(27), 5593–5599 (2002). [CrossRef] [PubMed]

24. G. Giovanelli, E. Palazzi, A. Petritoli, D. Bortoli, I. Kostadinov, F. Margelli, S. Pagnutti, M. Premuda, F. Ravegnani, and G. Trivellone, “Perspectives of 2D and 3D mapping of Atmospheric Pollutants over Urban Areas by means of airborne DOAS Spectrometer,” Ann. Geophys. 49(1), 133–142 (2006).

25. M. Premuda, S. Masieri, D. Bortoli, I. Kostadinov, A. Petritoli, and G. Giovanelli, “Evaluation of vessel emissions in a lagoon area with ground based Multi axis DOAS measurements,” Atmos. Environ. 45(29), 5212–5219 (2011), doi:, http://dx.doi.org/10.1016/j.atmosenv.2011.05.067. [CrossRef]

26. E. Palazzi, A. Petritoli, G. Giovanelli, I. Kostadinov, D. Bortoli, F. Ravegnani, and S. S. Sackey, “PROMSAR: A backward Monte Carlo spherical RTM for the analysis of DOAS remote sensing measurements,” Adv. Space Res. 36(5), 1007–1014 (2005). [CrossRef]

27. T. Wagner, J. P. Burrows, T. Deutschmann, B. Dix, C. von Friedeburg, U. Frieß, F. Hendrick, K.-P. Heue, H. Irie, H. Iwabuchi, Y. Kanaya, J. Keller, C. A. McLinden, H. Oetjen, E. Palazzi, A. Petritoli, U. Platt, O. Postylyakov, J. Pukite, A. Richter, M. van Roozendael, A. Rozanov, V. Rozanov, R. Sinreich, S. Sanghavi, and F. Wittrock, “Comparison of box-air-mass-factors and radiances for Multiple-Axis Differential Optical Absorption Spectroscopy (MAX-DOAS) geometries calculated from different UV/visible radiative transfer models,” Atmos. Chem. Phys. 7(7), 1809–1833 (2007). [CrossRef]

28. E. Cupini, M. G. Borgia and M. Premuda, Il codice PREMAR per la simulazione Montecarlo del trasporto della radiazione nell’atmosfera. (RT/INN/97/5. 1997).

29. E. Cupini, G. Ferro, N. Ferrari, Monte Carlo Analysis of Radiative Transport in Oceanographic Lidar Measurements. (ENEA/RT/INN/2001/7, 2001).

30. A. Berk, L. S. Bernstein, and D. C. Robertson, “MODTRAN: a moderate resolution model for LOWTRAN7,” Geophysics Laboratory, Air Force Systems Command, United States Air Force, Hanscom AFB, MA (1989).

31. F. Kneizys, D. Robertson, L. W. Abreu, P. Acharya, G. P. Anderson, L. S. Rothman, J. H. Chetwynd, J. E. A. Selby, E. P. Shettle, W. O. Gallery, A. Berk, S. A. Clough, and L. S. Bernstein, “The MODTRAN 2/3 Report and LOWTRAN7 MODEL,” Phillips Laboratory, Geophysics Directorate, Hanscom AFB, MA (1996).

32. P. Acharya, S. M. Adler-Golden, G. P. Anderson, A. Berk, L. S. Bernstein, J. H. Chetwynd, et al., “Modtran version 3.7/4.0 user’s manual,” Air Force Research Laboratory, Space Vehicles Directorate, Air Force MATERIEL Command, Hanscom AFB, MA (1998).

33. A. Berk, G. P. Anderson, P. K. Acharya, J. H. Chetwynd, L. Bernstein, and E. P. Shettle, “MODTRAN4 user’s manual,” Air Force Research Laboratory, Space Vehicles Directorate, Air Force Materiel Command, Hanscom AFB, MA (1999).

34. D. Nowak, L. Vuilleumier, C. N. Long, and A. Ohmura, “Solar irradiance computations compared with observations at the Baseline Surface Radiation Network Payerne site,” J. Geophys. Res. 113(D14), D14206 (2008), doi:. [CrossRef]

35. G. C. Pomraning, Radiation Hydrodynamics, Los Alamos National Laboratory Radiation Hydrodynamics Short Course Attendees, (1982), http://www.osti.gov/energycitations/servlets/purl/656708-zO5SF0/webviewable/656708.pdf

36. R. L. Liboff, Kinetic Theory: Classical, Quantum and Relativistic Descriptions. (Prentice Hall, 1989).

37. S. Chandrasekhar, Radiative Transfer, (Dover Publications, Inc., 1960).

38. S. Chatigny, M. Morin, D. Asselin, Y. Painchaud, and P. Beaudry, “Hybrid Monte Carlo for photon transport through optically thick scattering media,” Appl. Opt. 38(28), 6075–6086 (1999). [CrossRef] [PubMed]

39. C. P. Robert and G. Casella, Monte Carlo Statistical Methods (Springer, 2004).

40. A. De Matteis and S. Pagnutti, “Long-range Correlation-analysis of the Wichmann-Hill random number generator,” Stat. Comput. 3(2), 67–70 (1993). [CrossRef]

41. A. De Matteis and S. Pagnutti, “Controlling correlations in parallel Monte Carlo,” Parallel Comput. 21(1), 73–84 (1995). [CrossRef]

42. M. Premuda, “Montecarlo simulation of Radiative transfer in Atmospheric environments for problems arising from remote sensing measurements,” in Theory and Applications of Monte Carlo Simulations, S. Mordechai, ed. (INTECH, 2011).

43. M. H. Kalos and P. A. Whitlock, Monte Carlo Methods. Volume I: Basics (John Wiley & Sons, 1986).

44. M. Iqbal, An Introduction to Solar Radiation (Academic Press, 1983).

45. M. Premuda, S. Masieri, D. Bortoli, F. Margelli, F. Ravegnani, A. Petritoli, I. Kostadinov, G. Giovanelli, and E. Cupini, “A Monte Carlo simulation of radiative transfer in the atmosphere applied to ToTaL-DOAS,” Proc. SPIE 7475, 74751A, 74751A-12 (2009). [CrossRef]

46. K. Y. Kondratyev, Radiation in the Atmosphere (Academic Press, 1969).

47. F. X. Kneizys, S. A. Clough, E. P. Shettle, L. S. Rothmann and R. W. Fenn, “Linear Absorption and Scattering of Laser Beams,” Air Force Geophysical Laboratory, AFGL-TR-84 0265 (1984).

48. F. X. Kneizys, E. P. Shettle, W. O. Gallery, J. H. Chetwynd, L. W. Abreu, J. E. A. Selby, S. A. Clough and R. W. Fenn, “Atmospheric Transmittance/Radiance: Computer Code LOWTRAN6,” AFGL-TR-83–0187, (1983).

49. E. Palazzi, “Retrieval of trace gases vertical profile in the lower atmosphere combining Differential Optical Absorption Spectroscopy with radiative transfer models,” PhD thesis (Bologna University, Italy, 2008), http://amsdottorato.cib.unibo.it/983/1/Tesi_Palazzi_Elisa.pdf

50. E. Cupini, PREMAR-2 A Monte Carlo code for radiative transport simulation in atmospheric environments. (RT/INN/98/20. 1998).

51. A. Sarkissian, H. K. Roscoe, and D. J. Fish, “Ozone Measurement by Zenith-Sky Spectrometers: an Evaluation of Errors in Air-Mass Factors calculated by Radiative Transfer Models,” J. Quant. Spectrosc. Radiat. Transf. 54(3), 471–480 (1995). [CrossRef]

52. J. C. Chai, H. S. Lee, and S. V. Patankar, “Ray Effect and False Scattering in Discrete Ordinates Method,” Numer. Heat Transfer, Part B 24(4), 373- 389 (1993), http://www.informaworld.com/10.1080/10407799308955899

53. E. Papandrea, E. Arnone, E. Castelli, and M. Premuda, “Sounding the Upper Troposphere-Lower Stratosphere, with Satellite measurements,” Proc. “ESA Living Planet Symposium,” Bergen, Norway, 28 June- 2 July 2010 (ESA-SP-686, December 2010).

		Single scattering	Total scattering	Ground reflected	Total scattering plus ground
		radiance	radiance	radiance	reflected radiance
α^a = 180°	MODTRAN	2.12E-08	2.54E-07	1.78E-07	4.32E-07
ψ^b = 0°	MOCRA	2.12E-08	2.61E-07	1.78E-07	4.39E-07
θ^c = 30°	Diff. (%)	0.0	2.7	0.0	1.6
α^a = 180°	MODTRAN	9.09E-09	1.36E-07	7.75E-08	2.14E-07
ψ^b = 0°	MOCRA	9.07E-09	1.37E-07	7.74E-08	2.15E-07
θ^c = 60°	Diff. (%)	−0.2	0.7	−0.1	0.4
α^a = 150°	MODTRAN	2.03E-08	2.81E-07	1.58E-07	4.39E-07
ψ^b = 90°	MOCRA	2.03E-08	2.88E-07	1.57E-07	4.45E-07
θ^c = 30°	Diff. (%)	0.0	2.5	−0.6	−1.4
α^a = 150°	MODTRAN	1.04E-08	1.57E-07	6.87E-08	2.26E-07
ψ^b = 90°	MOCRA	1.04E-08	1.59E-07	6.85E-08	2.28E-07
θ^c = 60°	Diff. (%)	0.0	1.3	−0.3	0.9

		Single scattering	Total scattering	Ground reflected	Total scattering plus ground
		radiance	radiance	radiance	reflected radiance
α = 180°	MODTRAN	1.16E-08	6.85E-08	3.95E-07	4.64E-07
ψ = 0°	MOCRA	1.17E-08	6.77E-08	3.90E-07	4.57E-07
θ = 30°	Diff. (%)	0.9	−1.2	−1.3	−1.5
α = 180°	MODTRAN	6.75E-09	3.93E-08	1.97E-07	2.37E-07
ψ = 0°	MOCRA	6.78E-09	3.84E-08	1.95E-07	2.33E-07
θ = 60°	Diff. (%)	0.4	−2.7	−2.3	−1.7
α = 150°	MODTRAN	1.16E-08	7.68E-08	3.87E-07	4.64E-07
ψ = 90°	MOCRA	1.17E-08	7.53E-08	3.82E-07	4.57E-07
θ = 30°	Diff. (%)	0.9	−2.0	−1.3	−1.5
α = 150°	MODTRAN	1.16E-08	7.68E-08	3.87E-07	4.64E-07
ψ = 90°	MOCRA	1.17E-08	7.53E-08	3.82E-07	4.57E-07
θ = 60°	Diff. (%)	0.9	−2.0	−1.3	−1.5

		MOCRA		MODTRAN 4		Differences
	Solar source	Mean		Mean	Mean	Standard
λ (nm)	(watt/cm²·cm⁻¹)	Radiance		Radiance	Difference (%)	deviation (%)
310	5.1080E-07		6.272E-10	6.179E-10	2.4		4.7
435	3.2927E-06		6.131E-08	6.171E-08	1.3		2.6
550	6.1134E-06		1.648E-07	1.668E-07	−0.2		1.3
700	7.2381E-06		2.446E-07	2.502E-07	−1.8		0.5

		Single scattering	Total scattering	Ground reflected	Total scattering plus ground
		radiance	radiance	radiance	reflected radiance
α^a = 180°	MODTRAN	2.12E-08	2.54E-07	1.78E-07	4.32E-07
ψ^b = 0°	MOCRA	2.12E-08	2.61E-07	1.78E-07	4.39E-07
θ^c = 30°	Diff. (%)	0.0	2.7	0.0	1.6
α^a = 180°	MODTRAN	9.09E-09	1.36E-07	7.75E-08	2.14E-07
ψ^b = 0°	MOCRA	9.07E-09	1.37E-07	7.74E-08	2.15E-07
θ^c = 60°	Diff. (%)	−0.2	0.7	−0.1	0.4
α^a = 150°	MODTRAN	2.03E-08	2.81E-07	1.58E-07	4.39E-07
ψ^b = 90°	MOCRA	2.03E-08	2.88E-07	1.57E-07	4.45E-07
θ^c = 30°	Diff. (%)	0.0	2.5	−0.6	−1.4
α^a = 150°	MODTRAN	1.04E-08	1.57E-07	6.87E-08	2.26E-07
ψ^b = 90°	MOCRA	1.04E-08	1.59E-07	6.85E-08	2.28E-07
θ^c = 60°	Diff. (%)	0.0	1.3	−0.3	0.9

		Single scattering	Total scattering	Ground reflected	Total scattering plus ground
		radiance	radiance	radiance	reflected radiance
α = 180°	MODTRAN	1.16E-08	6.85E-08	3.95E-07	4.64E-07
ψ = 0°	MOCRA	1.17E-08	6.77E-08	3.90E-07	4.57E-07
θ = 30°	Diff. (%)	0.9	−1.2	−1.3	−1.5
α = 180°	MODTRAN	6.75E-09	3.93E-08	1.97E-07	2.37E-07
ψ = 0°	MOCRA	6.78E-09	3.84E-08	1.95E-07	2.33E-07
θ = 60°	Diff. (%)	0.4	−2.7	−2.3	−1.7
α = 150°	MODTRAN	1.16E-08	7.68E-08	3.87E-07	4.64E-07
ψ = 90°	MOCRA	1.17E-08	7.53E-08	3.82E-07	4.57E-07
θ = 30°	Diff. (%)	0.9	−2.0	−1.3	−1.5
α = 150°	MODTRAN	1.16E-08	7.68E-08	3.87E-07	4.64E-07
ψ = 90°	MOCRA	1.17E-08	7.53E-08	3.82E-07	4.57E-07
θ = 60°	Diff. (%)	0.9	−2.0	−1.3	−1.5

MOCRA: a Monte Carlo code for the simulation of radiative transfer in the atmosphere

Abstract

1. Introduction

2. Monte Carlo solution of the radiative transfer equation

3. Description of MOCRA

3.1 Radiation source and detectors

3.2 Geometrical and physical description of atmospheric environments

3.3 Path tracking and variance reducing techniques

3.4 Perturbative calculations

3.5 Computed quantities

4. Results, applications and comparison exercises

5. Summary and conclusions

Acknowledgments

References and links

Cited By

Figures (5)

Tables (3)

Equations (16)

Optics Express