SciELO - Scientific Electronic Library Online

 
vol.40 número1A robust algorithm for binarization of objects índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados

Articulo

Indicadores

  • No hay articulos citadosCitado por SciELO

Links relacionados

  • En proceso de indezaciónCitado por Google
  • No hay articulos similaresSimilares en SciELO
  • En proceso de indezaciónSimilares en Google

Bookmark


Latin American applied research

versión impresa ISSN 0327-0793

Lat. Am. appl. res. vol.40 no.1 Bahía Blanca ene. 2010

 

ARTICLES

An improved scheme for solving atmospheric radiative transfer problems with the spectral nodal method

M.P. De Abreu

Department of Computational Modeling, State University of Rio de Janeiro, P. O. Box 97282, 28610-974,
Nova Friburgo RJ, Brazil. deabreu@iprj.uerj.br

Abstract — In this article, we report on recent advances in a spectral nodal method for the numerical solution of multislab atmospheric radiation problems. Here, we derive a set of periodic relations for the coefficients of the bidirectional functions of our method, and we use these periodic relations to improve a recently developed computational scheme for solving a set of multislab atmospheric radiation problems with an arbitrary number and type of optically stationary layers. We present numerical results for a set of prototype problems that show the effect of stratospheric ozone depletion on the amount of ultraviolet-B radiation that reaches the Earth's surface.

Keywords — Multislab atmospheric radiation; Discrete ordinates; Periodic relations; Stationary layers; Ultraviolet-B radiation.

I. INTRODUCTION

In recent years, we have developed a spectral nodal method for numerically solving discrete ordinates radiative transfer problems in a stratified plane-parallel (multislab) grey planetary atmosphere (de Abreu, 2004a, b). In order to cast the discretized equations of our method in a suitable form for computer implementation and solution of practical problems in an efficient manner, we have successfully derived bidirectional functions for an arbitrary layer of a multislab model atmosphere (de Abreu, 2005a, b, 2006a). These bidirectional functions give us the emerging diffuse intensities at the bottom and top edges of a layer in terms of the incident diffuse intensities at top and bottom, and in terms of corresponding source terms. Moreover, they allow us to: i) replace corresponding atmospheric layers in multislab radiative transfer computations without loss of numerical accuracy; and ii) make an efficient use of our spectral nodal method for solving multislab atmospheric radiation problems with an arbitrary number and type (internal and/or boundary) of optically stationary layers, i.e. those layers in a multislab model atmosphere whose optical properties do not change from one problem to another in the set (de Abreu, 2006b).

In this article, we report on recent advances in our spectral nodal method. In Section II, we describe the class of radiative transfer problems that we are concerned with, and we outline our spectral nodal method. In Section III, we derive in full a set of periodic relations for the coefficients of the bidirectional functions of our spectral nodal method. In Section IV, we give a detailed account of an improved scheme for solving a set of multislab atmospheric radiation problems with optically stationary layers. In Section V, we present numerical results for a set of prototype problems in the ultraviolet-B (UV-B) range, and in Section VI we conclude this article with a discussion.

II. PROBLEM FORMULATION AND SPECTRAL NODAL METHOD

A. Problem Formulation

Let us assume that the multislab model atmosphere referred to at the onset consists of R contiguous, disjoint, and uniform layers so that layer 1 is the uppermost layer and layer R is the lowermost one. We assume further that the diffuse component of the intensity of the radiation field in the multislab domain can be accurately described by the discrete ordinates radiative heat transfer equations (Chandrasekhar, 1950; Siewert, 2000; Stamnes et al., 1988; Thomas and Stamnes, 1999)

(1)

We should note that the direct component of the radiation field is rather straightforward to model and compute (Thomas and Stamnes, 1999). In the set of Eqs. (1), N is the order of the angular quadrature, m from 1 to N/2 holds for the downward directions of radiation propagation, and m from N/2 +1 to N holds for the upward ones; t is the optical variable, and tr-1 and t r denote, respectively, the optical depths of the upper and lower boundaries of the rth layer; the symbol µm indicates the cosine of the angle defined by the discrete direction m and the downwardly oriented t axis, and the quantity ?n is the angular weight for direction n.

In this article, we use Gauss-Legendre quadrature sets (Lewis and Miller Jr., 1993), and we have ordered the discrete directions so that µm-1 < µm , m from 2 to N/2, and µm+N/2 = -µm m from 1 to N/2. Accordingly, ? n+N/2 = ?n, n from 1 to N/2, for the angular weights; is the single-scattering albedo in layer r; (2?+1)ß?,r is the ?th-order component of the Legendre expansion of the scattering phase function everywhere in layer r, which has been truncated after (Lr+1) terms; P? denotes the ?th-degree Legendre polynomial; the quantity µ0 is the cosine of the solar zenith angle, and 2p;µ0F0 is the top-of-atmosphere (t0) solar flux. As usual, we assume continuity conditions for the diffuse intensity at layer interfaces, and we use the discrete boundary conditions

(2)

B. Spectral Nodal Method

The spectral nodal method is a hybrid analytic/numerical method for solving the discrete ordinates Eqs. (1-2) without optical truncation error. That is to say, the numerical values generated on a digital computer for layer-average and layer-edge diffuse intensities are exactly those that would be obtained from the analytic solution of the discrete ordinates Eqs. (1-2), apart from computational finite arithmetic considerations and regardless of the optical properties of the layers. The method has been validated through a number of standard atmospheric radiation problems posed by the Radiation Commission of the International Association of Meteorology and Atmospheric Physics. More details can be found in de Abreu (2004a,b). From a practical viewpoint, the spectral nodal method is a discretization method that generates two sets of optically discretized equations without optical truncation error. The first set consists of the standard radiative heat balance equations

(3)

where ?tr tr-1 is the optical thickness of layer r, and are the rth layer-edge diffuse intensities in direction m at top and bottom, respectively; is rth layer-average diffuse intensity in direction m, and

(4)

is the rth layer-average first-collided source for direction m. The second set consists of the nonstandard auxiliary equations

(5)

A detailed description of the direct numerical schemes employed for determining the coefficients ?r,m,p and source terms gr,m can be found in de Abreu (2004b). Equations (3) and (5) can be reformulated as follows. We substitute Eqs. (5) into Eqs. (3) for corresponding layers. After some algebra, we obtain the set of N equations

(6)

where

(7)

and the first-collided source term

(8)

can be inferred from Eqs. (3) and (5). Upon half-range splitting (µm > 0 and µm < 0) of Eqs. (6), and after some algebra, we can write the bidirectional functions

(9)

and

(10)

where

(11)

and dm,p is the Kronecker delta.

The bidirectional functions (9) and (10) give us the emerging diffuse intensities at bottom (, m from 1 to N/2) and top (, m from N/2+1 to N) edges of layer r in terms of the incident diffuse intensities at top (, p from 1 to N/2) and bottom (, p from N/2+1 to N), and in terms of corresponding first-collided source terms for all r from 1 to R. Since the spectral nodal Eqs. (3) and (5) are free from optical truncation error, and optical approximations were nowhere introduced along the derivation steps, the bidirectional functions (9) and (10) are also free from optical truncation error. In this sense, the bidirectional functions (9) and (10) can be thought of as "exact discrete ordinates replacements" for corresponding layers throughout the multislab model atmosphere. The system formed by the bidirectional functions (9) and (10) together with the boundary conditions (2) can then be solved for the diffuse intensities at all layer edges and discrete directions with the one-cell block inversion iterative scheme (de Abreu, 2004b).

III. PERIODIC RELATIONS

We have shown (de Abreu, 2005c) that the coefficients ?r,m,p in the auxiliary Eqs. (5) do satisfy the periodic relations

(12)

and

(13)

for all r. In this section, we prove that the same relations are actually satisfied for all r, m, and p by the coefficients Hr,m,p in the bidirectional functions (9) and (10). We begin by showing that the coefficients hr,m,p given by definition (7) do satisfy the periodic relations (12) and (13). Since definition (7) holds for all r, m, and p, we may write

(14)

Since µm+N/2 = -µm, m from 1 to N/2, we now substitute the periodic relation (12) into (14) for the last term on the right hand side of expression (14) and find

(15)

We split the quadrature formula on the right hand side of Eqs. (15) into two terms so that

(16)

We now use the parity relation P?(µ) = (-1) ? P?(-µ), -1 = µ = 1, for the Legendre ?n+N/2 = ? n, n from 1 to N/2, for the angular weights, and the periodic relations (12) and (13) for the coefficients ?r,m,p to get

(17)

or, equivalently,

(18)

which is the periodic relation (12) for hr,m,p for all r. Starting from definition (7), we may write

(19)

Following similar steps from Eqs. (14) through (18), we get

(20)

which is the periodic relation (13) for hr,m,p for all r. We are now in position to show that the coefficients Hr,m,p given by definition (11) do satisfy the periodic relations (12) and (13). From definition (11), we may write

(21)

For m and p from 1 to N/2, we have dm+N/2,p+N/2 = dm,p, µm = |µm+N/2|, and hr,m+N/2,p+N/2 = hr,m,p. So, we can reformulate expression (21) further as

(22)

This is just the periodic relation (12) for the coefficients Hr,m,p for all r. Now, from definition (11) we may write

(23)

For m varying from 1 to N/2 and p from N/2+1 to N, we have µm = |µm+N/2|, and hr,m+N/2,p+N/2 = hr,m,p. So, we can reformulate expression (23) further as

(24)

This is just the periodic relation (13) for the coefficients Hr,m,p for all r, and our proof is complete.

It is worth noting that the coefficients Hr,m,p in the bidirectional functions (9) are exactly those on the right hand side of the periodic relations (22) and (24), while the coefficients Hr,m,p in the bidirectional functions (10) are exactly those on the left hand side of the periodic relations (22) and (24). Once the coefficients Hr,m,p in the bidirectional functions (9) have been computed, the coefficients Hr,m,p in (10) become readily available from the periodic relations (22) and (24). Therefore, the periodic relations (12) and (13) for ?r,m,p together with their companion relations (18) and (20) for hr,m,p, and (22) and (24) for Hr,m,p can be used to prevent unnecessary (and expensive) computation of RN2/2 coefficients ?r,m,p in the auxiliary Eqs. (5), plus RN2/2 coefficients hr,m,p in expression (7), and yet RN2/2 coefficients Hr,m,p in the bidirectional functions (10).

IV. AN EFFICIENT SCHEME FOR PROBLEMS WITH STATIONARY LAYERS

In this section, we describe an efficient scheme for solving on a digital computer a set of multislab atmospheric radiation problems with an arbitrary number and type of optically stationary layers. In this scheme, we make use of the spectral nodal method of Section II and the periodic relations of Section III.We begin by noting that the coefficients ?r,m,p in the auxiliary Eqs. (5) depend upon the quadrature dataset (µm and ?m for all m) and the optical properties of layer r (de Abreu, 2004a, b). Therefore, the coefficients hr,m,p and Hr,m,p behave similarly, i.e. they are uniquely determined by the above quantities, viz definitions (7) and (11). The quantities , which appear in the bidirectional functions (9) and (10), depend in addition on the cosine of the solar zenith angle (µ0). Now, suppose that we are given a set of multislab atmospheric radiation problems with an arbitrary number and type of optically stationary layers, and that each problem in the set is to be solved with our spectral nodal method for fixed N. By definition, the optical properties of the stationary layers do not change from one problem to another in the set, and so the coefficients Hr,m,p for all stationary layers can be computed only once and used in the appropriate pair of bidirectional functions for solving all problems in the set with our spectral nodal method. We do not have to recompute the coefficients Hr,m,p for the stationary layers each time a new problem in the set is to be solved. If in addition the solar zenith angle does not vary along the problems in the set, then the quantities for the stationary layers do not change from one problem to another, and can also be computed once for use in all problems.

We have updated our spectral nodal code (de Abreu, 2006b) written in FORTRAN 77 for solving a set of multislab atmospheric radiation problems characterized by an arbitrary number and type of optically stationary layers. The up-to-date version takes into account the periodic relations of Section III and follows the general lines of the above paragraph. For all stationary layers, we compute and store only half of the coefficients ?r,m,p, hr,m,p, and Hr,m,p by virtue of the periodic relations of Section III. Then, we divide the set of problems into a number of subsets, each of which formed by problems having the same value of µm (these subsets are what we used to call µ0 clusters). For a µ0 cluster, we compute and store the quantities appearing in the corresponding bidirectional functions (9) and (10) for all stationary layers. With appropriate exponential terms, the bidirectional functions can now be used to replace the corresponding stationary layers in all problems in the current µ0 cluster, and we solve one at a time with the one-cell block inversion iterative scheme. We then take to another µ0 cluster, and we proceed similarly until all problems in the set have been solved. We remark that if the solar zenith angle does not vary along the problems in the set, then the (only) µ0 cluster coincide with the problem set.

V. NUMERICAL RESULTS

In this section, we present numerical results for a set of prototype problems that show the effect of stratospheric ozone depletion on the downward flux of UV-B radiation that reaches the Earth's surface. For then we consider a four-layer model atmosphere (R = 4) over a dark surface with a single stationary layer (r = 1) for the upper atmosphere (?t1 = 0.04, = 0.3), a single stratospheric layer (r = 2) that we allow to change, and two stationary layers (r = 3 and 4) for the upper troposphere (?t3 = 0.07, = 0.4) and for the planetary boundary layer (?t4 = 0.1, = 0.64). A greytone representation of this model atmosphere is depicted in Fig. 1.


Figure 1. A four-layer model atmosphere.

We should note that the optical thickness and the single scattering albedo of the upper troposphere were obtained with the help of a specific-purpose UV radiative transfer code developed at the Norwegian Institute for Air Research (Engelsen and Kylling, 2005). We assume that the aerosol in the stratosphere consists of spherical (Lorentz-Mie) particles with an asymmetry factor of 0.9 in the UV-B range. The Legendre scattering phase function data for the uppermost layer (L1 = 8), the upper troposphere (L3 = 8), and the planetary boundary layer (L4 = 82) were extracted from a work of Garcia and Siewert (1985). The Legendre scattering phase function data for the stratospheric layer (we set L2 = 100 for simplicity) can be computed from the Legendre expansion of the Henyey-Greenstein phase function model for the aerosol particles (Liou, 2002).

From the model atmosphere, we define and solve a set of sixteen problems for the downward surface UV-B flux. Here, each problem is defined by the solar zenith angle and by the optical properties of the stratospheric layer in the UV-B range, which are mostly influenced by the stratospheric column ozone amount. We solve these problems with our spectral nodal code and with a computer code based on conventional numerical schemes for comparison. This computer code for comparison is based on the unconditionally stable diamond difference (DD) discretization scheme and on iteration on the scattering source, which are widely used schemes in neutron and photon transport computations (Lewis and Miller Jr., 1993). We remark that the numerical results in this section were generated with an IBM-compatible PC (2.67 GHz-clock Pentium D processor, 1 Gbyte of RAM) running on a GNU/Linux platform.

The sixteen problems were divided into four µ0 clusters. Each µ0 cluster consists of four problems corresponding to distinct pairs (?t2,?2 ). The stratospheric optical thicknesses shown in the first column of Table 1 were estimated from the aerosol optical thickness (a fixed contribution) and from the ozone absorption cross-sections in the UV-B (an average over the Hartley band and the smoothest part of the Huggins bands) for the following stratospheric column ozone amounts (in Dobson units): 320 (a typical value for tropical regions), 288 (a 10% decrease), 256 (a 20% decrease), and 224 (a 30% decrease) (Bais et al., 2006; Orphal, 2003). Therefore, the optical thicknesses indicated in the first column of Table 1 (from top to bottom) represent increasing levels of stratospheric ozone depletion.

Table 1. Downward surface UV-B flux (Wm-2).

aSpectral nodal results; bDD results.

In Table 1, we present numerical results generated by both the spectral nodal code and the DD code using N = 100 for the total (direct + diffuse) downward surface UV-B flux for each problem in the set. The numerical results in Table 1 show the expected trend toward increasing values of downward surface UV-B flux for increasing levels of stratospheric ozone depletion. It is worth noting that the numerical results generated by the spectral nodal code agree to within ±1 in the last figures given with corresponding DD results, which were generated on a very fine mesh grid of 50 optical cells per layer.

In reference to computer resources, we provide in Table 2 total computer memory location and execution time for solving the sixteen problems using the DD code (first row of results) and the spectral nodal code (second row of results). It is apparent that the spectral nodal code is more efficient than the conventional DD code to generate the desired quantities with comparable accuracy (note the remarkable reductions in both computer memory and execution time in the third row of results). This owes essentially to a couple of reasons. Firstly, the numerical results generated by the spectral nodal code for downward UV-B fluxes are free from optical truncation error. In order for the DD code to generate corresponding numerical results with comparable accuracy, it must use a very fine mesh grid in the computations. This leads the DD code to use up a considerable amount of computer resources (memory and time) in the computations. Secondly, the iteration on the scattering source scheme used in the DD code converges slowly for problems with at least one layer with a high value of single scattering albedo (Lewis and Miller Jr., 1993), and this is the case of the boundary and stratospheric layers for the sixteen problems. We remark that the execution time displayed in Table 2 for the DD code was obtained with the one-parameter Chebyshev extrapolation (Nakamura, 1986) option enabled to speed up the scattering source iterations.

Table 2. Computer memory and execution (CPU) time.

We close this section by noting that the variability in stratospheric ozone is not the dominant factor affecting the amount of UV-B radiation that reaches the Earth's surface, but rather the solar zenith angle; this dominance is indicated by the rows of Table 1. When the solar zenith angle is small (the sun is high, and µ0 approaches unity), the light optical path through the atmosphere is small, and so attenuation is lowered. Other factors affecting surface UV-B fluxes include seasonal variations in the Sun-Earth distance, altitude, aerosol distribution, surface conditions, and air pollution (Bais et al., 2006).

VI. DISCUSSION

We have reported on a number of improvements in a spectral nodal method for solving multislab atmospheric radiation problems. The extension of the periodic relations originally developed for the coefficients of the auxiliary equations to the coefficients of the bidirectional functions have led to significant savings in both computer memory and execution time, as half of the RN2 coefficients ?r,m,p in the auxiliary Eqs. (5), half of the RN2 coefficients hr,m,p in expression (7), and the coefficients Hr,m,p in the bidirectional functions (10) are no longer computed and stored. We remark that the computation of RN2 coefficients ?r,m,p involves the solution of N systems of linear algebraic equations in N unknowns each on a layer-by-layer basis. We have incorporated these periodic relations into a scheme developed by the author (de Abreu, 2006b) for solving a set of multislab atmospheric radiation problems with stationary layers. As a result, the up-to-date version of our spectral nodal code is able to solve these problems with increased efficiency. To illustrate these new features, we have considered a set of practical problems in the agenda of contemporary atmospheric science: the effect of stratospheric ozone depletion on the downward surface flux of UV-B radiation. To check the numerical results for accuracy, we have used a computer code based on the unconditionally stable diamond difference discretization scheme and on iteration on the scattering source, which are well-established, widely used schemes in radiation transport computations for their relative simplicity, level of discrete ordinates adequacy, and ease of computer implementation.

In closing, we should note that the incorporation of both specularly and diffusely (Lambertian) reflecting boundary conditions (Thomas and Stamnes, 1999) for the planetary surface into the current version of our spectral nodal method is underway. We intend to give a comprehensive account of the major changes in the bidirectional functions due to the presence of a reflecting planetary surface together with numerical results in a forthcoming article.

AKNOWLEDGMENTS

The work by the author was supported in part by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, Brazil).
The author is indebted to two anonymous referees for their helpful comments and constructive criticisms.

REFERENCES

1. Bais, A.F., D. Lubin, A. Arola, G. Bernhard, M. Blumthaler, N. Chubarova, C. Erlick, H.P. Gies, N. Krotkov, K. Lantz, B. Mayer, R.L. McKenzie, R.D. Piacentini, G. Seckmeyer, J.R. Slusser and C.S. Zerefos, "Surface ultraviolet radiation: past, present, and future", Report of the 2006 Assessment of Ozone Depletion, Chapter 7, United Nations Environment Programme, Ozone Secretariat, Nairobi, Kenya (2006).        [ Links ]

2. Chandrasekhar, S., Radiative Transfer, Clarendon Press, Oxford (1950).        [ Links ]

3. de Abreu, M.P., "A two-component method for solving multislab problems in radiative transfer", J. Quant. Spectrosc. Radiat. Transfer, 85, 311-336 (2004a).        [ Links ]

4. de Abreu, M.P., "Mixed singular-regular boundary conditions in multislab radiation transport", J. Comput. Phys., 197, 167-185 (2004b).        [ Links ]

5. de Abreu, M.P., "Layer-edge conditions for multislab atmospheric radiative transfer problems", J. Quant. Spectrosc. Radiat. Transfer, 94, 137-149 (2005a).        [ Links ]

6. de Abreu, M.P., "On efficiently solving a class of atmospheric radiative transfer problems using discrete ordinates models for bidirectional functions: the uppermost layer", Atmos. Sci. Lett., 6, 84-89 (2005b).        [ Links ]

7. de Abreu, M.P., "On increasing the efficiency of a discrete ordinates radiative transfer method with periodic relations", Therm. Eng., 4, 181-189 (2005c).        [ Links ]

8. de Abreu, M.P., "On equivalent diffuse conditions for an internal layer in multislab atmospheric radiation", J. Quant. Spectrosc. Radiat. Transfer, 102, 519-529 (2006a).        [ Links ]

9. de Abreu, M.P., "On the solution of multislab atmospheric radiation problems with an arbitrary number and type of stationary layers", Atmos. Sci. Lett., 7, 47-51 (2006b).        [ Links ]

10. Engelsen, O. and A. Kylling, "Fast simulation tool for ultraviolet radiation at the Earth's surface", Opt. Eng., 44, 041012 (2005).        [ Links ]

11. Garcia, R.D.M. and C.E. Siewert, "Benchmark results in radiative transfer", Transp. Theory Stat. Phys., 14, 437-483 (1985).        [ Links ]

12. Lewis, E.E. and W.F. Miller Jr., Computational Methods of Neutron Transport, 2nd ed., American Nuclear Society, LaGrange Park, IL (1993).         [ Links ]

13. Liou, K.-N., An Introduction to Atmospheric Radiation, 2nd ed., Academic Press, New York (2002).        [ Links ]

14. Nakamura, S., Computational Methods in Engineering and Science, Krieger Publishing Company, Malabar, FL (1986).        [ Links ]

15. Orphal, J., "A critical review of the absorption cross-sections of O3 and NO2 in the ultraviolet and visible", J. Photochem. Photobiol., A, 157, 185-209 (2003).        [ Links ]

16. Siewert, C.E., "A concise and accurate solution to Chandrasekhar's basic problem in radiative transfer", J. Quant. Spectrosc. Radiat. Transfer, 64, 109-130 (2000).        [ Links ]

17. Stamnes, K., 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, 2502-2509 (1988).        [ Links ]

18. Thomas, G.E. and K. Stamnes, Radiative Transfer in the Atmosphere and Ocean, Cambridge University Press, New York (1999).        [ Links ]

Received: September 29, 2008.
Accepted: December 7, 2008.
Recommended by Subject Editor Walter Ambrosini.