{'fraction_non_alphanumeric': 0.02780211901537481, 'fraction_numerical': 0.009049629570620568, 'mean_word_length': 4.962948888466654, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 11, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monte Carlo Simulation of Radiative Transfer in Atmospheric Environments for Problems Arising from Remote Sensing Measurements

The way in which solar radiation distributes itself in the atmosphere and on the ground is well known. It is beyond the scope of this book and the reader can refer to more specific references (Kondratyev, 1969; Goody & Yung, 1995; Liou, 1998) for more detail. Solar radiation, essentially in the visible-ultraviolet frequency range, and infrared radiation, emitted by the terrestrial surface, are the prevailing energy sources for general atmospheric circulation. They are thus particularly important for meteorological and climatic studies. It would therefore be of great interest, for instance, to be able to calculate the influence of the presence of ozone and trace gases, water vapour and clouds, and various aerosols on radiative transfer and global thermal energy in the atmosphere or in particular regions of it. These considerations naturally lead to the analysis of radiative transfer in the terrestrial atmosphere. This can be done using an atmospheric radiative transfer model (RTM) which also includes the possibility of single and multiple scattering events. Numerical and analytical methods can be used to solve a radiative transfer equation (Stamnes et al., 1988; Lenoble, 1977; Fouquart et al., 1980). A Monte Carlo approach is particularly suitable when multiple scattering significantly affects the results or where marked anisotropy of scattering and complex geometrical configurations are involved. The interest in such problems has increased through recently developed techniques related to remote sensing observations (satellite-based, ground-based or airborne) of the Earth's surface and atmosphere, involving the use of spectral radiation dispersion systems (mainly radiometers, spectrometers and interferometers) and active systems (principally RADAR, LIDAR and SODAR) Among these various techniques, DOAS (Differential Optical Absorption Spectroscopy) and LIDAR (Light Detection And Ranging) investigations on the presence of particular atmospheric constituents or of atmospheric phenomena such as clouds, fog, rain, etc., are of special interest. With reference to surface remote sensing observations, for instance, the effects of atmospheric absorption and scattering constitute a noise element, which has to be evaluated by calculations. Simulation of both LIDAR and DOAS systems deals with radiation in the UV/visible spectral range. For simulation purposes, a LIDAR system can be schematized as a pulsed

laser "disk" source of monochromatic radiation and a "disk" receiver (the "disk" schematization will be explained later), respectively with small emitting and receiving angles, i.e., with a narrow field of view (FOV). By means of Monte Carlo simulations, it is possible to analyze the time and spatial distributions of the backscattering radiation due to various atmospheric components. The equation for a LIDAR backscattering signal as well as further considerations relating to Monte Carlo simulations has been published (Pace et al., 2003;Coletti & Fiocco, 1980). DOAS (Noxon, 1975;Platt et al., 1979;Platt & Perner, 1980;Roscoe et al., 1999) is an established remote sensing technique 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 (passive DOAS). For simulation purposes, DOAS systems can be represented by a disk detector of solar radiation, with a relatively small diameter and a narrow FOV. For our purposes, the atmosphere can be considered to consist of two different kinds of components: molecular (gases) and non-molecular (aerosol), which interact with radiation in different ways. For both components, a knowledge of their interaction coefficients and of the angular distributions of the scattered radiation (phase functions) is required. Appropriate average values of interaction coefficients for suitable altitude subdivisions have to be computed by taking into account the height variation of their density. Regarding molecular components, both continuum and line absorption are considered. The former accounts for reciprocal interactions between molecules of the same or other species whereas line molecular absorption is connected to rotational and rotational-vibrational transition frequencies, characteristic of each kind of molecule. Using theoretical models it is possible to find analytical expressions for line absorption coefficients as functions of such frequencies (Clough, et al., 1981;Kneizys et al., 1983a;Kneizys et al., 1984;Kondratyev, 1969). For radiation scattering by molecular components, the well-known expression of scattering coefficients and phase functions arising from Rayleigh theory can be adopted (Kondratyev, 1969). For non-molecular atmospheric components, scattering coefficients and phase functions can be derived from the Mie theory (Kondratyev, 1969;Lenoble, 1977) for particles whose dimensions are comparable to the radiation wavelength. Besides the exact Mie phase functions, several approximate formulae are available, the most suitable being the Henyey-Greenstein approximation (Kondratyev, 1969;Lenoble, 1977), based on a knowledge of the asymmetry factor g (average cosine of the scattering angle). As can be seen in Figure 1 (Tomasi & Paccagnella, 1986), the Mie phase functions may show marked anisotropy. In the simulation of radiation transport through the atmosphere the phenomenon of refraction has to be taken into account as a consequence of different refractive index values between contiguous geometrical shells. This phenomenon can be remarkable, for instance, in the case of vertical view detectors receiving solar radiation at solar zenith angles near to or greater than 90° (horizon) as can be seen in Figure 2, where a plot of a sun ray, obtained using the MOCRA (MOnte Carlo Radiance Analysis) code for a 93° solar zenith angle (Premuda et al., 2009) is shown. This phenomenon is taken into account according to Snell's refraction law sin sin with respect to the normal to the boundary surface, respectively. It should be noted that total reflection occurs for incident angles greater than a critical value.
Z-Lin EARTH Fig. 2. Plot of refraction of a sun ray obtained using the MOCRA code for a 93° solar zenith angle. The refracted ray (dark blue) hits the Earth's surface (brown), whereas the linear ray, without refraction (magenta), would reach the vertical above the observer.
At ground level, the albedo phenomenon usually has to be taken into account according to appropriate surface albedo coefficients. Towards this end, the Lambert reflection law is frequently adopted, which assumes a cosine distribution law for the reflected radiation with www.intechopen.com respect to the normal to the ground surface. The values of the albedo coefficient may significantly depend upon the radiation wavelength and on the surface characteristics (snow cover, foliage, bare soil, etc.). Bidirectional reflectance distribution functions (BRDF), particularly useful for analyzing ground reflectance properties, may also be considered, as they take into account the solar zenith angle, the observer line of sight and the azimuth between the solar and observer directions. A theoretical model of the BRDF compared to experimental observations is available (Walthall et al., 1985). A Monte Carlo simulation has also been used (Richtsmeier & Sundberg, 2009). Very accurate molecular and non-molecular models are adopted in the radiancetransmittance MODTRAN codes (Berk et al., 1989;Kneizys et al., 1996;Acharya et al., 1998;Berk et al., 1999), where physical phenomena involving radiation in the infrared-ultraviolet spectral range are analyzed in detail. They are used widely by the remote-sensing community to model spectral absorption, transmission, emission and scattering characteristics of the atmosphere. Through a knowledge of interaction coefficients and phase functions, it is possible to set up sequences of discrete and continuous probability distribution functions (p.d.f.) which allow Monte Carlo photon-history tracking. Monte Carlo simulations related to LIDAR and DOAS systems will be described below, together with proper variance reducing techniques. A variety of Monte Carlo applications have been described (Marseguerra & Zio, 2002).

General features of Monte Carlo radiative transfer simulations in atmospheric environments
An appropriate description of a numerical or statistical simulation of radiation transport in an atmospheric environment requires a theoretical equation which governs the transport phenomenon to be stated and the physical and geometrical properties of the environment defined. To this end, in what follows, a simplified form of the integral transport equation will be given together with a possible environmental representation suitable for the simulation process, which will be subsequently described in its general fundamental lines.

Integral radiative transfer equation and Monte Carlo simulation
To analytically describe, in its photonic representation, radiation transport in atmospheric environments for time independent problems, the following integral form of the Radiative where k is the total extinction coefficient, sum of the scattering s k and the absorption a k coefficients. Q is the source density for emission and scattering defined as  (2) represents an internal radiation source. As can be seen, the specific intensity of radiation in direction Ω at point r is determined by the sum of two contributions: the first is the direct source contribution, due to the source specific intensity of radiation starting at point r s and the second accounts for emission and scattering of the beam radiation coming from each element ' ds of the path from s r to r along the Ω direction. Both contributions at point r are exponentially attenuated as a consequence of absorption and scattering collisions. It should be noted that in solar radiance analysis the first contribution accounts for the irradiance and the second for single and multiple scattering. The integral form of the transport equation is suitably solved by means of numerical or Monte Carlo simulations (Marseguerra & Zio, 2002). It can be obtained directly from the more general integro-differential radiative transfer equation (Spiga et al., 1992;Premuda & Palestini, 1982;Premuda, 1994). The attenuating exponentials are the same used to evaluate path lengths between collisions in Monte Carlo simulations. The length relevant to exponential attenuation of both contributions is the optical depth, defined as follows: In the standard "forward" Monte Carlo simulation, histories of photons emitted from the radiation source are followed until they hit the detector or disappear from the system as a consequence of leakage or absorption. Problems in this simulation procedure inevitably arise when small detectors with narrow FOVs are considered in the presence of an external broad source, due to the very low probability of a photon reaching the detector. To overcome this kind of problem, a "backward" Monte Carlo simulation is usually performed, taking into account that the linearity of the transport Boltzmann equation and the reciprocity relationship between the Green function G(P → P d ) and its adjoint G*(P d → P) make it possible to write for the flux Φ(P d ): : in which the source S(P) is always evaluated at the source point P (a detailed analysis of the Monte Carlo simulation of the adjoint transport Boltzmann equation is available (De Matteis & Simonini, 1978a;De Matteis & Simonini, 1978b)). In the case of solar sources, each photon history is therefore traced through the atmosphere as being generated by the detector according to its line of sight direction, taking into account the phase function for the angle which leads from a collision point P c to the source point P on the external atmospheric boundary along the solar radiation direction. A backward simulation description for radiation transport in the atmosphere, where polarization effects are also taken into account, is available (Collins et al., 1972). In atmospheric radiative transfer simulations, the backward scheme is even more suitable, because in photon interactions with atmospheric constituents energy changes do not occur. A forward Monte Carlo radiative transfer simulation for photon tracing in three-dimensional cloudy atmospheres is foreseen in the Mystic code (Mayer & Kylling, 2000). Topography and an inhomogeneous surface albedo are considered. www.intechopen.com

The atmospheric environment and photon history tracking
For a radiation of assigned wavelength, the photonic interpretation allows the application to radiative transfer of simulation techniques usually adopted for particle transport in assigned materials. Starting from the source, trajectories of individual photons are followed, according to the physical and geometrical properties of the environment of interest and taking into account the discrete and continuous probability distribution functions belonging to possible events for the photon. A synthetic description of the essential characteristics of the atmospheric environment required for the simulation process and of the standard "forward" history tracking is given below.

a) The atmospheric environment
Regarding the physical properties, as mentioned previously, reaction coefficients and phase functions for both molecular and aerosol particles are required. A knowledge of these parameters, together with refractive indices, which characterize each kind of atmosphere, can be profitably used to build a library data set, then used as the source of the fundamental physical and geometrical values needed for the Monte Carlo simulation. In analyzing radiative transfer in atmospheric environments, it can be crucial to verify the effects on radiation transport of perturbations in some atmospheric constituents (e.g. small variations in ozone or carbon dioxide concentrations). In Monte Carlo simulations it is possible to evaluate these effects by considering simultaneously reference and perturbed environments. This can be done by tracing the photon histories in the reference environment, taking into account the perturbations using appropriate weighting functions. This simulation tool makes it possible to evaluate small effects which could be masked by statistical errors when using separate calculations. This is possible, for instance, in the MOCRA code where several perturbed environments can be considered simultaneously. In radiation transport simulations in atmospheric environments, great advantage can be gained by proper representation of general 3D spherical multi-region geometries. The atmosphere can, for instance, be subdivided into cones, spherical shells and azimuthal half planes, obtaining very detailed geometrical descriptions of the various surface and atmospheric regions (Cupini et al., 2005). This allows one to take into account various ground altitudes and atmospheric profiles for each region, thus allowing for soil orography and latitudinal and longitudinal variations in the composition of the atmosphere. This is also possible with MOCRA and has been used to simulate the presence of an obstacle in horizontal passive DOAS measurements (Premuda et al., 2009). When several kinds of atmosphere have to be simultaneously taken into account in the same calculation, as may occur, for instance, in general 3-D geometries, the data belonging to each of them must obviously be available in the library. The appropriate build-up and treatment of the physical and geometrical library data can make it possible, for instance, to foresee effects due to the injection of special aerosols (e.g., fumes) in a given atmospheric region with already assigned aerosols. In the following Km are the units of length and, consequently, Km -1 are reaction coefficient units. As already pointed out, reaction coefficients, given by the product of the density and the proper particle reaction "cross section", will depend on the particle density behaviour along the vertical z-axis starting from the ground altitude, so that average values are to be computed in Monte Carlo simulations. To this end, the total assumed atmospheric height www.intechopen.com range (z 0 , z max ) can be subdivided into an assigned number NL of conveniently chosen geometrical layers. With reference to molecular particles, if one assumes an exponential behaviour of the density spatial distribution within each layer, depending upon the pressure and temperature profiles, one obtains, for the j-th layer: The average density value j within the layer will be given by: i.e.: If, for numerical reasons, this relationship cannot be applied, the value j can be assumed to be equal to the arithmetic mean of j and j+1 . Similar relationships will hold, consequently, for the average molecular scattering and absorption Ma j k coefficients. It can be observed that, although a unique value is usually given for the scattering coefficient, corresponding to that of air considered as a fictitious molecule with a molecular weight M Air = 28.964, average absorption coefficients must be assigned for each molecular species. The average coefficient Ma j k has to be interpreted as the sum over the average coefficients belonging to each of them. Analogous considerations can be made for the average aerosol coefficients As j k and Aa j k for the same layers, and for the refraction indices. It can be seen that very different kinds of aerosols can be found along the z-axis, from the boundary layer to the high stratosphere, so that marked discontinuities can occur in the average reaction coefficients between contiguous layers. For the j-th vertical layer, the following total coefficients can be considered: According to eq. (3), the optical depth of the layer will be: jj j =k h. This optical depth characterizes the layer in the atmospheric environment and gives a measure of the transparency of the layer to the incident radiation. For the sake of simplicity, in what follows the brackets will be omitted. In the photon diffusion process, an essential role is played by molecular Rayleigh scattering, which affects the distance travelled and the motion direction following a collision. A theoretical formula for the scattering coefficient Ms k which takes into account the dependence on wavelength , refraction index n and molecular number density N is given by Kondratyev (Kondratyev, 1969) together with phase function and refraction index expressions. A refraction index formula for standard air, which takes into account the wavenumber dependence, is available (Edlén, 1966). In the LOWTRAN-MODTRAN codes the following formula is adopted (Kneizys et al., 1983b): 1.70 10 c= ⋅ . A molecular scattering formula, which takes into account the depolarization coefficient δ , is given by (Kneizys et al., 1984) It is stressed that the Rayleigh scattering coefficient rapidly increases as the radiation wavelength decreases. For the depolarization coefficient, a value of 0.0279 is assumed, for instance, in the MODTRAN code. In this code, the following formula is used for the phase function () θ p (Kneizys et al., 1983b): With reference to the aerosol reaction coefficients and phase functions, as mentioned above, the Mie theory can be applied, which holds for spherical homogeneous particles with radii r comparable to the wavelength λ of the incident radiation. The Mie theory is much more complex than Rayleigh theory, from both analytical and numerical points of view, and will not be described here. It will only be observed that, for a particle of radius r , the size parameter 2 r/ = and complex index of refraction m=n iq − are quantities of interest in the theory. If a distribution

( )
Nr is available (see, for instance, Shettle, 1989 for possible distribution functions), average reaction coefficients over the distribution can be obtained. Regarding the phase function, approximate equations have been proposed. These include the well-known Henyey-Greenstein phase function, given by where g , the asymmetry factor, is the average cosine of the scattering angle, i.e.: In MODTRAN, special models for different kinds of rain droplets are also considered. The aerosol phase functions can be made available in the library by means of tables of values over assigned angles in the range (0° -180°). For the molecular phase functions, the formula given by equation (11) can be used directly during the simulation. Attention must be paid in sampling scattering angles: for molecular scattering, the scattering angle cosine can be sampled using equiprobable cosine tables; for aerosol scattering, on the other hand, such tables are, as a rule, not possible to set up, due to the marked anisotropy in their phase functions (Premuda, 1994).

b) Photon history tracking
Once the physical and geometrical properties of the environments required for the simulation process have been established, the history tracking can be carried out as usual, taking into account the peculiarities of the problem to be solved. In this regard, it can be observed that in the efficiency of the simulation aimed at the evaluation of a given quantity a crucial role is played by the estimator chosen, which will assume a statistical value i at the end of i-th history. If N histories have been processed, its average value : will give the estimate of the quantity of interest. As for the statistical error, this can be estimated by means of the standard deviation for both the estimator, ε σ and the estimate, , through the well-known variance equations 2 V= and 2 V= : A quantity used to evaluate the efficiency of the simulation process adopted is given by the so-called Figure of Merit (FOM): where t is a measure of the computation time. As an example of the various possible choices of the estimator in the simulation process, one can consider the problem of finding the absorption by a given molecule, for instance ozone, of the radiation at an assigned wavelength in a layer. Two kinds of estimators can be adopted: the first one could be derived from the probability of the radiation being absorbed in a collision within the layer (collision estimator); the second by multiplying the distance travelled by the photon within the layer by the absorption coefficient of the molecule (distance estimator). This latter estimator can be conveniently applied when very few collisions occur within a frequently crossed layer.
To take into account in the simulation process the analytical information available on the physical events which may occur and to develop tools for enhancing the contribution of the photon history to the final result, a statistical "weight" w is associated with the travelling photon, taking an initial value 0 w assigned at the beginning of the history. In the solar radiance analysis, for instance, 0 w could be assumed to be equal to the radiation intensity 0 I of the solar source. Regarding the photon source, the spatial and angular distributions of the emitted photons will make it possible to select for each history the starting parameters needed for the subsequent realization of the chain of possible events undergone by the photon in a given reference system. For this purpose, the importance must be stressed of the www.intechopen.com availability of a proper generator of random numbers uniformly distributed over the interval (0, 1), usually an arithmetic congruence, which, as much as possible, avoid internal correlations among the generated numbers (Knuth, 1981;De Matteis & Pagnutti, 1988). It is well known that discrete ( 1 ... From this equation, for instance, a value τ of the optical path travelled by the photon is easily obtained from x e − attenuating exponential function, normalized to 1 over the interval ( ) 0,∞ . When a direct solution is not available or too expensive, special techniques (such as the so-called "rejection technique") can be used to obtain values of the variable x obeying the p.d.f of interest. In any case, discrete cumulative probability distribution tables can be obtained by means of analytical or numerical integration for assigned values of the independent x variable: A value of x can be obtained from eq. (17), by means of proper linear interpolations for the chosen random number r. This is, for instance, the case of the Mie phase function, for which, as previously said, tables of values are assumed to be available for the diffusion angle or, more conveniently, its cosine. Taking into account, as already emphasized, the possible marked anisotropy of this distribution function, tabulating points are assumed to be accurately chosen to avoid improper choice of the scattered photon direction. Regarding the Rayleigh phase function, the same approach can be adopted. In this case, as already indicated, a given number (preferably a power of 2, e.g. 32) of cosine intervals can be derived, by means of interpolations, from the cumulative probability aerosol or by a molecule is made. If a collision with an aerosol particle occurs and various kinds of aerosols exist within the layer, each with its proper phase function, with probabilities:

A2s A2s
As k p= k , …, the specific kind of aerosol involved will be selected. On the basis of the proper diffusion cosine distribution function, a new motion direction will be chosen and a new path started. If the distance d to be travelled is greater than the distance from the starting point and the boundary of the layer along the motion direction and a leakage phenomenon from the geometrical system does not occur, the optical path exhausted within 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 1 n and 2 n belonging to the two contiguous layers, the new refracted motion direction will be computed. It must be taken into account that a total reflection will occur if 1 2 n> n and, moreover, the angle of the incident direction on the boundary between the two layers is greater than the L θ angle, being For termination of the history tracking, various criteria can be adopted in addition to that from the leakage of the photon from the geometrical system. Among them, the following can be considered. 1) Assigning as input parameter a minimum weight fraction, min w , the history will end if 0 min w <w w . The minimum weight fraction has to be carefully chosen to avoid expensive calculations giving non-essential contributions to the estimator. A value of 10 -5 could be, for instance, adopted in many cases. The sum over the histories of the unprocessed weights will give a measure of the bias introduced. "russian roulette" game is played: with a probability 0 r w p= w that the photon will survive and the initial weight 0 w=w will be restored; with a probability 1 r p − that the photon will be killed and the photon weight lost. At the end of all the histories processed, the sum of the total weight gained should be statistically equal to the total weight lost. To avoid large jumps in the weight, which can affect the variance, 0.5 cut w= could be a reasonable value to be assigned. In MOCRA, all three options are available.
According to the quantities to be estimated through the simulation process, analytical contributions can conveniently be computed during the history tracking, so as to reduce as www.intechopen.com much as possible the variance of the calculation. Regarding, for instance, the molecular absorption in a layer, an example is given by the distance estimator previously described. A further estimate frequently required is given by the leakage from the geometrical system. At each collision point the optical distance e to the external boundary along the direction of motion is computed. The escape weight e e w= w e − will be the analytical contribution of the collision point to the searched-for estimate. The total sum of the analytical contributions will give the estimator value for that history. In addition to analytical evaluations during the history tracking, variance reducing techniques can be devised which alter the natural sequence of the events undergone by the photon, locally or over the entire physical and geometrical system, giving rise to more efficient unbiased contributions to the required estimates. The most common is given by the well-known forced collision technique. If the direction of motion from a starting or a collision point does not cross a reflecting boundary, the leakage optical distance e is computed as previously described. The two possibilities for the photon escaping from the system with a probability , normalized to 1 over the interval ( ) 0, e , will then be adopted for a photon of weight c w=w . When this technique is applied, the first criterion among those previously described to end the history tracking has to be used. In this case, as a precaution, a maximum number of collisions per history can be assigned. It must be observed that the forced collision technique is an expensive one and it may cause, moreover, considerable changes in the travelling photon weight which can affect the variance. It should be used when the collisions within a system of small optical thickness are of special interest. In any case, the FOM should be taken into account. Further special variance-reducing techniques will be described within the context of the following LIDAR and DOAS simulation analysis. As previously highlighted, to analyze the effect of environmental perturbations on the final result due to small changes in physical parameters which do not affect the phase functions, a simulation process which takes into account, simultaneously, the unperturbed (reference) environment and the perturbed one can be adopted. As an example of application, the Air Mass Factor calculation (described in section 3.2) according to its theoretical definition can be considered, where the solar radiances with and without the absorption of a particular molecule are required. In its basic lines, the procedure can be derived from the so-called "importance sampling" technique: given an assigned p.d.f, ( ) f x , normalized to 1 over ( ) According to this device, a photon weight will be associated with the perturbed environment. The tracking of the photon history will be carried out in the reference system, taking into account the perturbed one by updating its weight according to proper weighting factors. As an example, if, in the reference system, along its line of flight a photon crosses a boundary, a weighting factor f w= e − , ' = − , has to be used, where and ' are the optical paths from the starting point to the boundary for the reference environment and for the perturbed one, respectively. If a collision occurs, the weighting factor will be can be assumed. Special care must be taken, in any case, over the treatment of the various events during the history tracking, such as, for instance, those connected with the forced collision technique. An estimator of the differential effects of the perturbation upon a given result can be directly set up by collecting, during the history tracking, the differences of interest between the perturbed and the unperturbed environments. The exponential dependence of the weighting fractions causes this instrument to be suitable for very small perturbations, difficult to evaluate using separate calculations. Obviously, several perturbed environments can be envisaged in the same calculation, each of them being associated with a statistical weight. It must be remarked that, according to this technique, the leading photon history is governed by the reference system, so that the end of its history will cause the end of the histories in the perturbed ones, too. If the min w criterion is assumed to end the history tracking, it is convenient, for each perturbed environment, to sum the total unprocessed weighting. The perturbation technique is foreseen in MOCRA (see below for applications) and in the successive versions of the PREMAR code. A special perturbative technique is also available (Rief, 1984). To handle statistically significant estimators and to avoid too many contributions being collected, which can cause a loss of precision, batches of assigned numbers of histories can be envisaged. For each batch, the average of the quantity of interest over the processed histories will give the estimator value belonging to that batch. The average value over the number of batches will give the searched estimate. As an example, 100 batches each of 1000 histories could be run. To make it possible to reach a desired precision in the simulation, a restart option over the batches is envisaged, so as to optimize the calculation time. To facilitate the comparison between two different calculations characterized by small differences in the physical or geometrical parameters, a strategy on the random number generation can be devised for both batches and histories inside a batch, so that each history will begin with its proper initial random number. Unaltered corresponding histories in the two calculations will give in this way the same results for the same quantity. This device is envisaged both in MOCRA and in PREMAR (Premuda et al., 2009;Cupini et al., 2001).

Simulation of LIDAR and DOAS observations
In general "Remote Sensing" techniques analyze kinds of interactions between a wave and the examined medium. For minor gases in atmosphere electromagnetic radiation is certainly www.intechopen.com the most adequate and its interactions with the atmosphere are known as "spectroscopic". Various spectroscopic techniques have been developed to measure concentrations of atmospheric pollutants. There are two main types of remote sensing systems: passive remote sensing and active remote sensing. Passive sensors detect natural radiation emitted or reflected by the object or surrounding area being observed. Reflected or scattered sunlight is the most common source of radiation measured by passive sensors. Examples of passive remote sensors include film photography, charge-coupled devices and radiometers. Active collection, on the other hand, emits energy in order to scan objects and areas whereupon a sensor then detects and measures the radiation that is reflected or backscattered from the target. RADAR is an example of active remote sensing where the time delay between emission and return is measured, establishing the location, height, speed and direction of an object. Remote sensing makes it possible to collect data on dangerous or inaccessible areas. Its applications include monitoring deforestation in areas such as the Amazon Basin, the effects of climate change on glaciers and on Arctic and Antarctic regions, and depth sounding of coastal and ocean depths. Remote sensing also replaces costly and slow data collection on the ground, ensuring in the process that areas or objects are not altered. Chemical remote sensing systems (radiometers, spectrometers, interferometers) are mounted on many satellites, allowing mapping of pollutants or of minor atmospheric components over large areas of the earth. Orbital platforms thus collect and transmit data from different parts of the electromagnetic spectrum, which in conjunction with aerial or ground-based sensing and analysis, provides researchers with enough information to monitor trends such as El Niño and other natural long-and short-term phenomena. Doppler radar is used, for instance, in enhanced meteorological collection such as wind speed and direction within weather systems. Other types of active collection include plasmas in the ionosphere. Interferometric synthetic aperture radar is used to produce precise digital elevation models of large scale terrain (see RADARSAT, TerraSAR-X, Magellan). Radiometers and photometers are the most common instruments in use, collecting reflected and emitted radiation over a wide range of frequencies. The most common are visible and infrared sensors, followed by microwave, gamma ray and rarely, ultraviolet. They may also be used to detect the emission spectra of various chemicals, providing data on chemical concentrations in the atmosphere. Among the various remote sensing techniques, LIDAR (Light Detection And Ranging) and DOAS (Differential Optical Absorption Spectroscopy), allowing the analysis of gas and aerosol concentrations in the atmosphere. The applications of DOAS remote sensing systems are numerous on the ground, and in both airborne and satellite configurations. LIDAR ground systems are used to detect and measure the concentration of various chemical compounds in the atmosphere, while airborne LIDAR can be used to analyze the ground characteristics, such as the heights of objects and vegetation, more accurately than with radar technology. Underwater LIDAR and LIDAR in a coupled air-ocean system allow the analysis of marine environments in the seas and oceans, for oil spills and, phytoplankton development. In water systems, fluorescence emission and Raman scattering can play a fundamental role. Below, the Monte Carlo simulation of LIDAR and DOAS observations is discussed, together with proper variance-reducing techniques. www.intechopen.com

LIDAR systems
A LIDAR system consists of a laser pulsed source at a given wavelength and a receiving telescope collecting the backscattered radiation. From the characteristics of the detected radiation (total intensity, time and spatial distributions) it is possible to obtain information about the nature and concentration of particles or about obstacles encountered along the path. Given an assigned reference system (Oxyz) with z-axis normal to the ground, for simulation purposes, the source S and the telescope T can be schematized as independent plane disks centered at given coordinate points in the geometric environment, each with its own diameter, axis and FOV. In Fig. 3 a graphical schematic 2-D representation is shown. Regarding the telescope, the disk representation can be adequate for atmospheric LIDAR systems, where collision points normally occur far from the telescope disk centre. More realistic descriptions of the telescope device are needed, for instance, in underwater LIDAR systems, as envisaged in PREMAR (Cupini et al., 2001). Fig. 3. Schematic disk representation of LIDAR system for simulation purposes. S is the source disk, T is the telescope disk, P1 and P2 are collision points, θ1 and θ2 are angles between flight direction and direction connecting collision points to the telescope centre.
In Monte Carlo simulation an outgoing ray from the vertex of a fictitious source cone, possibly located outside the atmospheric geometrical system, is uniformly generated, i.e., if α S is the angle which characterizes the source aperture angle, an angle cosine with respect to t h e c o n e a x i s i s u n i formly chosen in (cosα S, 1) with a corresponding azimuth angle uniformly chosen in (0, 2 ). The ray intersection with the source disk identifies the initial coordinates of the source photon, whose motion direction coincides with that of the ray. If the source aperture angle is equal to zero, the photon starting point is sampled uniformly on the disk with its motion direction orthogonal to it. For photons whose collision point falls into the telescope FOV the distance R between such a point and the centre of the telescope disk is evaluated (see Fig. 3). The expected collision contribution to the intensity T I of the radiation collected by the telescope can be written as: where as w is the product of the current photon weighting and the probability of scattering by aerosol particles, ms w the analogous quantity for a molecular component, ( ) a P Θ and ( ) m P Θ the corresponding phase functions for the angle Θ between the flight direction before scattering and the direction from the collision point to the centre of the telescope disk, normalized to 1 over the whole solid angle, is the optical distance between collision point and the centre of telescope disk. The solid angle element Ω is given by the expression where T A is the area of the telescope disk and Φ the angle between the telescope axis and the disk centre-collision point direction. To avoid infinite variances the calculation can be performed only for R values greater than a certain threshold. For this purpose, if small distances between the collision point and the telescope disk are of interest, as could happen, for instance, in underwater LIDAR systems, the following Ω formula can be used, where n r is the radius of the disk normal to the direction of the photon-to-telescope-disk centre: which is reduced to eq. (20) when /0 n rR → . In Fig. 4 a comparison between the approximate and the correct formulae is given. The aerosol A I and molecular M I contributions to the intensity T I of radiation in eq. (19) can be considered separately in the calculation. According to assigned time intervals, the time elapsed distribution from the emission at the source to detection at the telescope can be www.intechopen.com obtained. Along the axis of the telescope disk, starting from its surface plane, spatial intervals can be considered. Spatial distributions of the radiation intensities on the telescope can thus be obtained by projecting onto its axis the collision points within the FOV. Within this context, a possibly interesting quantity is the spatial distribution of the so-called backscattering ratio between the total and molecular radiation intensities. A comparison between the Monte Carlo results for the backscattering ratio obtained with PREMAR and experimental ones obtained at the ENEA centre at Brasimone lake (Bologna), located at an altitude of 0.91 Km, is given in Fig. 5. The scheme of the vertical Lidar system is given in Fig. 6 (not to scale), the disk source having a diameter d 1 = 2.8 cm with a FOV of α 1 = 0.12 mrad and a telescope disk diameter d 2 = 80 cm with a FOV = 0.35 mrad, with a distance of 80 cm between the two disk centers. A wavelength of 0.532 m was considered. The measurements were carried out over two periods: the first in September 1993, to verify the presence of volcanic aerosols following the eruption of Pinatubo, which occurred in June 1991; the second in February 1995. As can be seen, a satisfactory agreement is obtained when using the high volcanic profile foreseen by MODTRAN for the fall-winter season (on the left) and background stratospheric aerosol profile, again predicted by MODTRAN for fallwinter (on the right).
(a) (b) Fig. 5. Backscattering ratio R: comparison between PREMAR simulation and measurements carried out with LIDAR at the Brasimone ENEA centre in September 1993, after the eruption of Pinatubo, which occurred in June 1991, and in February 1995. The simulation was performed using: a) a high volcanic MODTRAN profile for fall-winter ; b) a background stratospheric aerosol MODTRAN profile for fall-winter (Cupini et al. 1997).
To reduce the variance of the calculation, besides the forced collision technique previously described, the local forced collision and the splitting techniques can be used when the contribution to the backscattered radiation of particular layers with small optical depths is of special interest.
According to the local forced collision technique, when a photon particle crosses a boundary of the layer of interest, the photon coordinates of the crossing point are memorized and a virtual collision is forced to occur in the layer with the proper statistical weight for the photon. The foreseen statistics related to this collision point are performed and the memorized coordinates of the photon redefined. A path is then chosen in the usual way, without further statistics if a collision occurs within the layer. In this latter case, the device can be repeated as in the case of the crossing of the layer boundary (see Fig. 7). The advantage of such a technique resides in the fact that efforts are concentrated only upon the layers of interest and, moreover, the photon travelling weight, unlike the general forced collision technique, remains unchanged. Regarding the splitting technique, when a photon crosses a boundary of the layer of interest with an assigned splitting index m, analogously to a local forced collision, the photon coordinates of the crossing point are memorized. From this point a number m of photons are generated, each with a weight w i = w/m (i=1,…,m). All the photons are processed within the layer, selecting their starting path according to the exponential distribution, ending their history, except for the last one, if they leave it. If the last photon emerges from the layer, its www.intechopen.com weight is multiplied by the splitting index m and its history tracking continues (see Fig. 8). If more adjacent layers with the same splitting index m are simultaneously considered, the processing of the generated photon can continue into the new layer with the same splitting index. Adjacent layers with different splitting indices can be considered, as foreseen in PREMAR but, for simplicity, are not described here.
Excessively high values of m should be avoided: a value of between 2 and 10 should be adequate. For special treatment of the splitting technique, see, for instance, Burn (Burn, 1995;Burn 1997). It is of interest to underline that the local forced collision and the splitting technique can be adopted within the same calculation for the same or different layers. Fig. 8. Splitting technique for a layer with splitting index m=4. The first photon undergoes a collision and its history ends when it hits the layer boundary. The second and the third photon would leave the layer without colliding but their histories end. Only the history of the last generated photon continues, as it emerges from the layer after a collision.

DOAS systems
DOAS is a 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 wavelength range. The molecular absorption is analyzed to obtain the concentration of the trace gases integrated along the optical path between the source and the receiver (Slant Column Density, SCD). After being introduced by Noxon (1975) during stratospheric trace gas studies at Fritz Peak observatory in Colorado, DOAS quickly became one of the most promising methods for determining the role of minor gases affecting ozone depletion in the lower stratosphere. In the late 1970s, Platt and Perner carried out DOAS measurements of CHOH, O 3 , and NO 2 in maritime and rural areas of Northern Europe (Platt et al., 1979;Platt & Perner, 1980). During the following decades DOAS techniques have been applied to urban pollution analysis and monitoring as well as in the climatic and environmental fields following two different methodological and instrumental approaches: 1) the active mode which can, by using an artificial source of radiation, perform measurements of atmospheric minor gas concentrations, integrated along the optical path between the lamp and the receiving telescope of the instrument (Perner et al., 1976;Evangelisti et al., 1995;Stutz & Platt, 1997a,b).
2) The passive mode, using diffuse solar radiation as the radiation source, carries out measurements of the examined gas column content along both the vertical (zenith-sky) and oblique (also called Off-Axis or MAX-DOAS (Wagner et al. 2007a)) directions. The zenith-sky configuration is especially used for observations in the lower stratosphere region for research mainly related to climate studies; the Off-Axis mode is mainly applied in environmental monitoring (Petritoli et al., 2002;Hönninger et al., 2004;Bortoli et al. 2005, Giovanelli et al., 2006. Particularly in the passive mode, both in vertical and oblique configurations, an Air Mass Factor (AMF) is calculated as a fundamental part of the interpretation of DOAS data in atmospheric observations using a Radiative Transmission Model (RTM) which also requires the use of a multiple-scattering configuration to describe the radiation passage through the atmosphere.
The AMF concept can be generalized to the recently developed "off-axis" configuration (Hönninger et al., 2004) and appropriate RTMs have to be used for data interpretation.
Another recent application of this type of remote sensing system is the so-called ToTaL-DOAS (Topographic Target Light scattering DOAS) which is a novel experimental procedure to retrieve trace gas concentrations present in the lower atmosphere. Scattered sunlight reflected from natural or artificial targets located at different distances are analyzed to retrieve the spatial distribution of the concentration of different trace gases such as NO 2 , SO 2 and others (Frins et al., 2006;Frins et al., 2008;Louban et al., 2008). In this case RTMs are required to compute the Equivalent Path Length (EPL), i.e. the distance of a fictitious radiation source giving the same signal collected by the detector (Premuda et al., 2009). For the simulation of radiometer or DOAS devices, a reference system O(x,y,z) can be assumed with its origin at the Earth's centre and the z-axis normal to ground surface. Direction cosines u s , v s , w s will define the solar source radiation direction. An external atmospheric boundary surface is considered, over which the starting points of the solar radiation are projected. Two kinds of detectors can be taken into account: the first being schematized as a point detector located on the ground, looking upward, and the second as a disk detector of a given radius, with its central point geometrical co-ordinates (x 0 ,y 0 ,z 0 ) and normal viewing direction (u 0 ,v 0 ,w 0 ). In the latter case, the FOV is assumed. The disk receiver can be reduced to a point with an assigned line-of-sight (u 0 ,v 0 ,w 0 ).
In the first case, if ( ) z I is the radiance intensity due to a solar ray reaching the detector from a height z, the total intensity T I can be obtained by integrating the contributions from all values of z from the ground to the upper boundary of the atmosphere: If the entire vertical range is subdivided into N intervals of thickness i z (i=1,…,N) it can equivalently be written as In the Monte Carlo simulation, a stratified sampling procedure can be performed, as in MOCRA, which consists of randomly choosing a set of z points uniformly distributed within each height interval i z , and estimating the corresponding ( ) Iz value. The average intensity value over the points belonging to the i-th interval times the vertical layer thickness i z gives the contribution of that interval to the total intensity. The single scattering radiance S I is calculated using a forward Monte Carlo simulation of the photon path from the solar source coordinate on the external geometrical boundary, computed on the basis of the z altitude of each selected point and of the assigned solar zenith angle, www.intechopen.com taking into account the refraction when the photon path crosses a boundary between geometrical regions with different refractive indices, until the z-axis is reached. The different arrival points on the z-axis will be utilized to compute the contribution to S I of the appropriate intervals. The average geometrical and optical slant and vertical paths can be evaluated. The multiple scattering radiance M I is calculated using a backward Monte Carlo simulation from the selected vertical points t o t h e s u n , a s p r e v i o u s l y d e s c r i b e d . T h e contribution of each interval to M I is computed and the total scattering radiance TSM I= I+ I evaluated. In the M I calculation, albedo and refraction phenomena can be envisaged, possibly excluding, as in MOCRA, refraction along the path from the last scattering point to the sun. Perliski and Soloman (Perliski & Solomon, 1993) carried out a backward Monte Carlo simulation in a spherical shell model atmosphere to calculate the air mass factor for a vertical upward-looking detector. In case of disk detectors with assigned diameter, FOV, position and orientation, let P d =(x d , y d , z d ) be the centre of the disk, D (≥ 0) its diameter and (u d , v d , w d ) the direction cosines of the outgoing direction normal to the disk. Moreover, let α (> 0) be the field of view of the detector. By assuming D > 0, a fictitious cone with a semi-amplitude α and the disk as its base can thus be generated. The vertex point P v ≡ (x v , y v , z v ) (possibly located outside the atmospheric geometrical system as previously mentioned for the LIDAR source) can be considered as the virtual source point in a backward Monte Carlo simulation, with initial motion directions uniformly distributed within the cone. For each Ω direction, the crossing point with the detector disk will give the true starting point for the particle. Let Ω be the solid angle amplitude and ( ) I Ω the radiance intensity on the receiver from the Ω direction. The total intensity T I over all possible directions is assumed to be: Eq. (24) can be simulated through a standard backward Monte Carlo estimating procedure for T I . 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 is that of the detector axis. The PROMSAR (PROcessing of Multi -Scattered Atmospheric Radiation) code (Palazzi et al. 2005) was developed to perform backward Monte Carlo simulations of DOAS observations for vertical and off-axis looking detectors. In MOCRA, through a 3D multi-region geometry, optionally selectable by the user, the possibility of defining different topographic and atmospheric scenarios for a set of user-defined regions is included. This allows simulating, for instance, a horizontal point detector with a line of sight perpendicular to a vertical obstacle 100 m high and 50 m wide, at a distance of 1 Km from the detector. The solar direction cosines were obtained by means of the astronomical parameters. Fig. 9 shows a comparison of the simulation results with measurements obtained from the DOAS spectrometer TropoGAS (Tropospheric Gas Analyzer Spectrometer) working, in this case, in a horizontal-view configuration at a distance of 1 Km from a house (Premuda et al., 2009 Fig. 9. Passive-active measurements comparison at S. Pietro Capofiume (11.6226° E, 44.6538° N), Bologna, Italy. Measurements were performed almost simultaneously by the same instrument in active (Xenon lamp) and passive modes. The lamp was placed in front of a house at 1 Km from the spectrometer. A passive measurement towards the house was inserted in the measurement table. Only data with errors of less than 30% are plotted. The axis covers the time interval 7-19 (hours) with a linear scale.

Sun-lamp comparison
Data obtained from passive measurements by means of the simulation show a very good agreement with the simultaneous active measurements used as a reference for the comparison. As for the AMF calculation, let I and I * be the radiance (single or total scattering) detected by the receiver with and without the trace gas of interest, respectively. The AMF of the species is defined as (Sarkissian et al., 1995): where a is the vertical absorption optical depth of the molecular species of interest. In the case of weak absorption, an approximate formula for AMF is given by : where δ OPT is the intensity-weighted absorption optical path over all the collisions or reflections which contribute to the total intensity I * . If II * ≈ can be assumed, the calculation can be performed directly within the reference system containing the trace gas (Sarkissian et al., 1995). In MOCRA, when calculating OPT , perturbed radiance values are used if perturbative calculations are performed, whereas unperturbed radiance values are used if this is not the case. Due to the possibility of perturbative calculations, an estimation of the correct air mass factor can be carried out considering a perturbed environment without the absorption contribution of the molecular particle of interest. Fig. 10 shows plots of exact and approximate AMFs for single-scattering radiance, computed by MOCRA for O 3 (Fig. 10.a) and NO 2 (Fig. 10.b) at wavelength of 310 nm for a mid-latitude summer atmospheric environment with urban extinction. It can be seen that for ozone there is a significant difference between the exact and approximate value at higher solar zenith angles, due to its marked absorption at the given wavelength enhanced by the greater path length, whereas the values for NO 2 , which is a weak absorber, are almost identical to each other. To generate the gas profiles, starting from the gas SCDs via an inversion method, the "box air mass factors" which describe the sensitivity of the measurements as a function of atmospheric layer altitude have to be computed. They are defined as (Pukıte et al., 2006): where b n is the number density of the gas in the box, b h is the box height and b β is the absorption coefficient, VCD b (Vertical Column Density) is the integral of concentration along the vertical direction within the box. This is computed, for instance in the TRACY (Trace gas RAdiative Transfer Monte Carlo Y(I)implementation) series models (Pukıte et al.,2006;Kühl et al., 2008), simulating radiative transfer, using a backward Monte Carlo scheme.
A Monte Carlo evaluation of the box AMF can be performed as previously described for the AMF over the entire atmospheric environment, by considering the box with and without the trace gas of interest and assuming the optical box height as a in eq. (25) and eq. (26). In eq. (26) the absorption optical distances which effectively contribute to the box AMF will be those belonging to the box. In MOCRA several boxes can be considered simultaneously, each consisting of one or more vertical layers.
PROMSAR Box Air Mass Factors and radiances were validated through comparison with a series of state-of-the-art UV/visible RTMs (Wagner et al., 2007b).
To optimize the type of calculation which involve single-and multiple-scattering radiances being computed together with, where envisaged, the reflected radiance, simulation techniques can be adopted which force the path starting from the source to undergo the desired event. In the case of single-and multiple-scattering radiances, for the photon, born at the source with an initial weight 0 w=w , the probability c p of collision within the system is computed. According to this probability, a photon with a starting weight 0c w=w p is forced to collide, thus giving its contribution to the desired radiances. If the reflected contribution from the ground to the detector is not required, the simulation can stop and a new source photon can be chosen. Otherwise, if the initially selected motion direction crosses the ground, a new photon, with the same characteristics of the previous one, starts from the source with a weight ( ) , which will be reflected according to the envisaged albedo coefficient. In the case where only the reflected contribution is of interest, this last simulation is the only one performed. It should be noted that, if scattered and reflected contributions are both required, each history is split into two sub-histories for a given source-starting photon. Such a technique can efficiently be applied when there are low probabilities that the starting photon will collide or be reflected. In MOCRA the standard simulation procedures together with those described here are envisaged. It should be stressed that this technique is independent from the general, previously described forced collision technique, which can not be applied when the direction of motion crosses the ground. The forced collision technique can in any case be applied for collisions following the first one or reflection. Clearly, such a device could also be adopted for Lidar system simulations, such as when ground characteristics are to be examined using airborne or satellite Lidar systems.

Possible future developments
The simulation methods presented here are essentially devoted to the solution of RTE in the UV-visible spectral regions, where the internal source term ( ) Sr in Eq. (2) can be neglected. But remote sensing problems, such as MIPAS (Michelson Interferometer for Passive Atmospheric Sounding) data analysis, dealing with the infrared spectral region, requires a consideration of the presence of sources at the Earth's surface and within the atmospheric system, so that ( ) Sr cannot be neglected. In fact, the atmosphere presents a wide absorption window in the visible frequencies, whereas it is a strong absorber in the infrared, and, according to Kirchoff's law, the emission and absorption coefficients for each wavelength are the same. The principal problem to be solved is thus how to represent these two kinds of sources. To a first approximation, the external source at the Earth's surface should be set equal to Planck's blackbody function ( ) www.intechopen.com In general, the internal source and the scattering and absorption macroscopic cross-sections depend upon the description of atoms and molecules which constitute matter and their energetic status and there is no simple relationship between them. One frequently advanced hypothesis is that of Local Thermodynamic Equilibrium (LTE) (Pomraning, 1982). This assumes that medium properties are dominated by atomic collisions and at each time at each point the system is in thermodynamic equilibrium and the radiation field, even though very different from ( ) B ,T , does not influence thermodynamic equilibrium. If this hypothesis holds, at each point of the system and at each time t , only two variables (besides the system composition) are required to calculate the internal source S and the absorption and scattering cross-sections a and s . Thus, with the usual meaning of the symbols, On the basis of these considerations a further development of the work done within this framework could be the extension of MOCRA to the infrared spectral range introducing a source at the Earth's surface defined by Eq. (28) and an internal radiation source defined by Eq. (30) according to surface or air temperature and radiation frequency. Blackbody internal sources are foreseen in MODTRAN codes. More detailed studies should be devoted to molecular scattering simulation, taking into account, besides the Rayleigh scattering, Raman scattering and their mutual connections (Young, 1981). In PREMAR, Rayleigh and Raman scattering are considered separately for water environments.