1. Introduction

The Earth’s atmosphere, which is essential for life, behaves as a perfect multilayered turbid medium with absorbents and scatter elements. When sunlight strikes the atmosphere, it can be reflected, absorbed, scattered or/and transmitted. Each of these phenomena occur depending on the atmospheric composition. The stratosphere contains ozone as a major absorbent element. The intensive absorption peaks of radiation by ozone occur in wavelengths that range from 240 to 320 nm. Simulating the propagation of light through the atmosphere which is a complex geometry medium is not easy. This is because the atmosphere has an extensive composition and undergoes continual changes. Currently, there is still a problem with calculating the optical parameters of the atmosphere which are the absorption coefficient, scattering coefficient, anisotropy factor, and refractive index.

The way in which the propagation of light in earth’s atmosphere is studied is through the Radiative Transfer theory. Such theory appears with the studies of electromagnetic radiation in the interstellar medium¹ which triggers the radiative transfer equation (ETR). The ETR determines the propagation of light through the specific intensity². However, it has not been trivial finding a deterministic solution to the ETR.

There are several analytical attempts to get a solution, such as the theory of diffusion approximation³, the simplified theory of Kubelka-Munk⁴, and numerical attempts, which are seen in the case of the Monte Carlo method⁵. In addition, a distinct method that has been tried in order to obtain the ETR solution is the decomposition of Fourier series⁶. However, the study in which the diffusion approximation is used has its limitations and is only an approximation to the solution of the ETR. This is because it only offers a solution through the differential equation of heat which has user restrictions. Similarly, the Kubelka-Munk theory is another approximation of the diffusion approximation as it only applies to specific geometric conditions. On the other hand, the numerical Monte Carlo method is suitable when it comes to multiple scattering environments with complex geometric configurations⁷. Ultimately, the Fourier series depicts another attempt to solve the ETR. Although it is mathematically complicated, it only concludes that it is possible to apply Fourier series in the radiative transfer calculations. However, a study that provides a concise solution through the use of this method does not exist.

The advancement in computers’ management of data has allowed to use numerical methods for the finding of a solution to the ETR. Furthermore, the use of Monte Carlo’s numerical method to obtain a solution for the ETR has allowed designing optical diagnostic techniques that are noninvasive.

In this regard, the numerical Monte Carlo method becomes relevant, and when coupled with the rise of computing power, a diagnostic tool emerges. Monte Carlo is a powerful technique used for solving complex problems. In this method, computers simulate the behavior of real world systems based on a statistical sampling theory and probabilistic analysis of complex physical systems⁸. In the ETR absorption ( μ 𝑎 ) and scattering ( 𝜇 𝑠 ) coefficients, and the anisotropy factor (𝑔)⁹, appear because they are optical parameters used in the simulations. These coefficients depend on numerous factors including the properties of the medium and the wavelength of the incident light. In simulations that utilize numerical Monte Carlo method, it is possible to define the rules of photon propagation in the form of probabilistic density functions which describe the probability of absorption or scattering. This provides a flexible solution to this problem and allows deriving a methodology to decode the information associated with radiation distribution through a turbid media such as the atmosphere.

Thus, behind all the management of numbers and data through Monte Carlo the idea is to retrieve the absorption coefficient of ozone in the stratosphere through the trigonometric adjustment of the diffuse reflectance spectra. The spectral data is adjusted with the trigonometric series using the Fourier expansion to find dependence between the terms of the expansion of the Fourier series and the absorption coefficient of ozone in the stratosphere. As a result, a method that allows relating the parameters of the adjusted spectral curves of diffuse reflectance with the optical properties of the stratosphere emerges.

2. Atmospheric Model

The National Aeronautics and Space Administration¹⁰ defined in its rules the U.S. standard atmosphere as a steady-state representation of the first 100 km from the ground where a moderate solar activity period is assumed. Remarkably, beginning at 51 km atmospheric tables are identical to the rules of the US Standard Atmosphere¹¹. However, the earth’s atmosphere is generally described in terms of its layers. Generally, the atmosphere is considered to be one that extends over 560 km on the planet’s surface and is divided into four layers the troposphere, stratosphere, mesosphere and thermosphere, with the respective interfaces: tropopause, stratopause and mesopause. However, what allows for the complete understanding of the atmospheric model is to denote the important processes that occur within it, such as absorption and scattering. The atmosphere is composed of small particles, which are air molecules, mixed with large particles at low pressures called aerosols. Given the continually evolving composition of the atmosphere the propagation of light in it is erratic and is shown in Fig. 1. The scattering of light by air molecules is studied by Rayleigh¹² while the scattering of light in the presence of aerosols is studied by Mie¹³.

To perform a simulation of the diffuse reflectance of light in the atmosphere, the Monte Carlo for Multi-layered media (MCML) program was modified. This allows to model a photon’s pathway in a multilayer medium that contains a refractive index (n) which depends on the respective layer where the optical properties ( 𝜇 𝑎 , 𝜇 𝑠 ,𝑔) and thickness (d) of each layer are used as input parameters to simulation¹⁴. The multilayer medium used is the proposed simplified atmosphere model (Fig. 2). It consists of two layers: the troposphere and stratosphere which are 15 km and 50 km in thickness respectively.

Consequently, the selection of absorption and scattering elements to consider in this study is closely related with the selected spectrum. That is, the absorption of radiation is a selective process with the wavelength primarily due to aerosols, clouds, and atmospheric gas components. Instead, scattering

Figure 1 Atmosphere Light Propagation. 

Figure 2 Atmospheric Model. 

is a non-selective process which means, it affects all wavelengths. Also, the atmospheric constituents responsible for this process are gaseous molecules, ice crystals, water droplets, and, like in the absorption process, aerosols. Scattering is a process that preserves the total amount of energy, but the direction in which it spreads can be altered. Absorption is a process that removes the energy of the electromagnetic radiation and converts it into another form of energy. And the extinction or attenuation is the sum of the scattering and absorption that occurs, which represents the effect on the radiation as it passes through a turbid medium.

The outermost layer of the atmosphere analyzed in the proposed model is the stratosphere. One of the main constituents of the stratosphere is ozone, which is considered the primary absorbent agent uniformly distributed therein. In this regard, the greatest pondered weight in the ozone is assumed as main absorbent. Thus, the absorption coefficient 𝜇 𝑎 in the stratosphere is calculated with Eq. (1);

(1)

That is, 90% absorption due to stratospheric ozone and the rest due to water vapor, nitrogen dioxide, and oxygen in equal parts.

The second layer discussed in the proposed model input is the troposphere. In the troposphere, water vapor, carbon dioxide, nitrogen dioxide, and oxygen are considered as absorbent agents. As a result, the absorption coefficient 𝜇 𝑎 in the troposphere is calculated according to Eq. (2) as shown below:

In other words, the absorption due to each of the elements in the troposphere is considered in equivalent parts.

In this situation, the main absorbents considered in the troposphere and stratosphere is¹⁵: water vapor (H 2 O), carbon dioxide (CO 2 ), oxygen (O 2 ), ozone (O 3 ) and dioxide nitrogen (NO 2 ).

The input data for the simulations with the modified MCML of absorption in the troposphere and stratosphere of each component were taken from the publication of Spectral Atlas which is of free access for the public on the website of the Institute of Chemistry, Max-Planck Mainz, in Germany¹⁶. Some of the data files were adjusted to the range of wavelengths of the study and others mathematically adjusted via Mat lab in order to obtain the applicable mathematical relationship, as shown in Table 1.

With regards to the input data of the simulations with the modified MCML, in terms of scattering, it is studied based on the size of the particle and it affects all wavelengths. The phenomenon of scattering or attenuation of radiation occurs when the constituents forming the atmosphere scatter the radiation passing through varying their initial propagation direction. According to the relationship between the size of the wavelength of the incident radiation and the media particles, it is possible to divide this scattering phenomenon¹⁷ in:

Rayleigh scattering, this phenomenon occurs when the wavelength of the incident radiation is much larger than the size of the scattering particles.

Mie Scattering, this phenomenon occurs when the scattering particles are of the same order of size as the wavelength of the incident radiation.

Geometrical optics, it occurs when the wavelength of the incident solar radiation is much smaller than the size of obstacles in its propagation.

Therefore, scattering, in this study, is calculated according to the approach based on the combination of Mie and Rayleigh¹⁸⁾ shown in Eqs. 3 and 4, where A and B represent the weighted fraction of relaxation in each layer.

Coming up next, the proposed model of the atmosphere as entry for the simulations (Fig. 2) is presented, In terms of the refractive index (n) in the atmosphere¹⁹ and the non-planar parallel geometry of such, it is observed that it strongly depends on the altitude, pressure, temperature and wavelength of the radiation measured. In turn, in emptiness 1 substitutes for the refractive index and as it approximates to the Earth’s surface 1.00029. The factor or parameter of anisotropy, which measures the degree of anisotropy of the scattering, is defined as the average of the scattering angle, which is obtained by the phase function, and its value varies in between -1 and 1²⁰. For particles with a radius in the range of 3-30 microns, 𝑔 values are in the range of 0.8 to 0.9. For aerosol particles whose radius’ are typically ∼0.1 𝜇m, 𝑔 values are less than 0.5∼0.7²¹. Therefore, the anisotropy factor is considered constant throughout the studied spectrum with the value of 0.6 in the stratosphere and 0.8 in the troposphere. Predicting the distribution of light inside any turbid media requires detailed information about the optical properties of the medium²². Distinguishably, the entry model

Table 1 Absorbents, data & math relations. 

presented in (Fig. 2) summarizes the optical behavior of a stable atmosphere.

3. MCML Simulation

It has been established that Monte Carlo is a relevant and flexible technique for simulating the propagation of light in turbid media, and that the optical study of materials and other components provide information of the surface morphology and internal structure due to the interaction of light with matter. This approach is useful because make it possible to propose a diagnosis of the medium in a noninvasive way.

In the simulation process through the modified MCML, the input is formed by the beam of photons (photon number N), the refractive index of the medium (𝑛), the absorption and scattering coefficients, the anisotropy factor (𝑔) and the distance (𝑑) or size layer to be analyzed. On the other hand, the output is formed by the diffuse reflectance.

The simulation is based on the aleatory pathway of the photons through the medium, which are selected statistically according to each absorption and/or scattering event. After propagating a large number of photons, the results are usually closer to reality. To simulate the diffuse reflectance, the Monte Carlo algorithm multilayers of Oregon Medical Laser Center, which is software of free access⁵ was modified. The modified algorithm requires the following input parameters: the absorption coefficient 𝜇 𝑎 , scattering coefficient 𝜇 𝑠 , anisotropy factor 𝑔, refractive index 𝑛, and the thickness 𝑑 of each of the layers in each of the simulations performed. The input parameters used in the simulations of the atmosphere are shown in Fig. 2. The simulations through the modified MCML allow obtaining the spectral curves of diffuse reflectance versus the wavelength.

4. Results

In order to retrieve the absorption coefficient of the ozone in the stratosphere the ozone concentration is varied from 1% to 10 % and the rest of the concentrations remain fixed. Therefore, changes in the spectral curves are attributed to a unique physical parameter. In the Fig. 3, a family of spectral curves of diffuse reflectance which were obtained from the simulation is shown as the concentration is varied by 2.

Figure 3 Curves of Diffuse Reflectance versusWavelength varying the Ozone Concentration. 

Evidently, it is difficult to see what happens to the diffuse reflectance when the ozone concentration increases. Thus, numerical parameters for these curves are determined by applying adjustments through trigonometric series. For each of the simulations, the trigonometric adjustment is performed using Fourier series grade 2 as shown in Eq. (5).

(5)

In addition, it is assumed that the spectral curves are periodic. Hence, the natural frequency from the definition of the range where the curve is periodic can be determined and as a result comply with the conditions of the Fourier series. The Eq. (6) shows the calculation of the natural frequency used.

The adjustment of the data with the least squares nonlinear method using Eq. (5) is shown in Fig. 4. The adjusted spectral curves of diffuse reflectance in the Fig. 4 show that the residual behavior is lower and more stable in comparison with the curves in Fig. 3.

Furthermore, the natural frequency is calculated using Eq. (6). It makes possible to obtain the Fourier coefficients that are 𝑎 0 , 𝑎 𝑛 and 𝑏 𝑛 , which adjust the simulated data.

Table 2 Fourier Coefficients versus % concentration of O3 

Figure 4 Adjusted Curves of Diffuse Reflectance versus Wavelength. 

Figure 5 Absorption Coefficient O3 versus Fourier Coefficient a1. 

That is, the diffuse reflectance increases or decreases depending on such parameters, making it possible to establish a relationship between the change in the concentration of ozone and the Fourier coefficients. Table 2 shows the parameterization constants for each of the simulated reflectance curves.

An analysis of the behavior of each constant of parameterization of the Fourier series was performed by varying the concentrations of O 3 . As a result, a chaotic behavior for most constant was obtained, except for constant 𝑎 1 , whose behavior can be observed in Fig. 5.

Thus, a relationship between the concentration of O 3 and the parameterization constant 𝑎 1 is established (Eq. 7). This, in turns, allows retrieving the concentration of O 3 from the following polynomial function.

(7)

5. Discussion

The input parameters in the proposed atmosphere model were defined with a polynomial function grade 5 for the Ozone. These input parameters were used as initial input for the simulation through the modified MCML. The simulations through the modified MCML allow obtaining the spectral curves of diffuse reflectance versus the wavelength as shown in Fig. 3, with a linear concentration of ozone changes. These results or diffuse reflectance curves obtained from the simulations were analyzed by Fourier series, which is a nonlinear fit. The spectral curves in Fig. 3 were parameterized by the trigonometric function, in order to obtain the coefficient 𝑎 1 which is related to the experimental data and associated with different concentrations of ozone in the stratosphere.

In this regard, to analyze the behavior of the coefficient 𝑎 1 while the ozone concentration varies, a polynomial grade 5 is obtained, which is similar to the polynomial ozone proposed as input in the atmospheric model. It should be noted that the polynomial of grade 5 governs the amount of ozone in the stratosphere, while it is impacted by the component 𝑎 1 of the Fourier series. However, despite being able to obtain the starting point for simulations, it is even more important to recognize the self consistency that is demonstrated when the Fourier coefficient 𝑎 1 varies in a similar way as the polynomial proposed in the input parameters for the Ozone absorption. That is, 𝑎 1 orderly varies with the reflectance changes allowing retrieving the absorption coefficient of stratospheric ozone. It should be noted that beyond obtaining a relationship (Eq. 7) that applies only to the model of the atmosphere proposed, a method that retrieves the optical parameters is obtained due to the sensitivity of the Fourier coefficient 𝑎 1 because of the changes in the absorption coefficient of stratospheric ozone.

6. Conclusions

The study presents a retrieving process of an optical parameter in the stratosphere through the simulated spectral curves of diffuse reflectance, as a new approach. This allows concluding that it is possible to obtain a relationship between the coefficient 𝑎 1 of the expansion of Fourier series and the absorption coefficient of the ozone in the stratosphere through trigonometric adjustment as an effective method.