Photographic observation and optical simulation of a pollen corona display in Japan Souichiro Hioki1,2 and Hironobu Iwabuchi1,* 1
Center of Atmospheric and Oceanic Studies, Graduate School of Science, Tohoku University, 6-3 Aramaki-aza-Aoba, Sendai 980-8578, Japan 2
Currently at Department of Atmospheric Sciences, Texas A&M University, MS 3150 TAMU, College Station, Texas 77843-3150, USA *Corresponding author: [email protected] Received 2 July 2014; revised 13 August 2014; accepted 14 August 2014; posted 15 August 2014 (Doc. ID 215247); published 2 October 2014
Brightness and chromaticity profiles were extracted from a vivid solar corona image taken with a digital camera in Sendai, Japan, to compare with a radiative transfer simulation applying Lorenz–Mie theory and single-scattering approximation. The comparison revealed suspended particles having a narrow particle size distribution peaking at radius 14.5 μm. Presumably, pollen of an indigenous coniferous tree, the cryptomeria (Cryptomeria japonica), is responsible for the corona display. The extracted brightness and chromaticity profiles are reproduced well by assuming the presence of a water soluble aerosol and dust in addition to the pollen. We find that photographic analysis of corona displays, similar to that used to measure cloud particle size, is applicable to estimating pollen particle size distribution and column number density. © 2014 Optical Society of America OCIS codes: (010.1290) Atmospheric optics; (010.1110) Aerosols; (010.1310) Atmospheric scattering; (010.1690) Color. http://dx.doi.org/10.1364/AO.54.000B12
1. Introduction
Coronas are iridescent concentric rings around the Sun or the Moon. This atmospheric optical display consists of multiple color-separated rings with apparent radii of a few degrees and is present when particles with radii of a few micrometers are uniformly suspended. As the radii of suspended particles can be estimated from the apparent size of a corona using diffraction theory [1], the estimated particle size is often reported with corona displays due to cloud particles and pollens. Observational approaches to coronas widely utilize photographic analysis in estimating particle radii [2–4]. Among them, Sassen (1991) [2] presents the largest number of samples. Some observed cloud 1559-128X/15/040B12-10$15.00/0 © 2015 Optical Society of America B12
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
coronas have noncircular shapes [3,5], implying effects of cloud heterogeneity or phase transition. Noncircular coronas can be caused by pollen. In 1992, a Finnish amateur astronomers’ network reported vertical elliptical coronas [6]. A simulation of elliptical coronas that uses the scattering theory of spheroids showed that nonspherical pollens can induce noncircular or partly brightened coronas, due to aerodynamic preferences of falling particle orientation [7]. The production of irregular coronal displays from particle orientation was confirmed by experimental simulation [8]. The shape and the size of pollen corona rings depend on the pollen species. Coronas caused by birch, pine, and juniper pollens have been reported [9]. Numerical simulations of pollen corona displays are limited because of difficulties in the calculation of scattering properties of nonspherical particles and radiative transfer with oriented particles.
However, cloud corona displays have become increasingly accessible over the past decade. As a result of increased computing power, scattering calculations for spherical particles by using Lorenz–Mie theory (Mie theory, below) are now possible even on personal computers. Laven showed that fast Mie theory can simulate rainbows, coronas, and glories [10]. Recent simulations of cloud corona displays realistically include background Rayleigh scattering, multiple scattering by clouds, the curvature effect of the Earth, and the effect of ozone absorption [11,12]. This is a major improvement from earlier works such as [13], which calculated scattering patterns by water cloud particles only. Pollen often causes serious allergic symptoms, and is, thus, a public health concern. In Japan, mature coniferous forests of cypress and cryptomeria (or Japanese cedar, or Sugi; Cryptomeria japonica), have been increasing since the 1960s as a result of governmental afforestation policies and the decline of Japan’s forestry industry in later years [14]. These trees cause allergy symptoms among about a third of the population of Japan. Cryptomeria pollen, mainly released from March to April of each year, has been well documented regarding size, allergens, and preferential weather conditions for release. Cryptomeria pollen is nearly spherical, except for a small bump on its top, with a radius ranging from 12 to 16 μm [15]. The pollen count in Sendai, Japan reportedly increases when the temperature increases and humidity concurrently decreases, especially after rain [16]. Despite atmospheric pollen having a great impact on human activities, there have been only limited optical analyses of pollen corona displays. To the best of the authors’ knowledge, the only photographic analysis of pollen coronas in Japan is that of Takahashi et al. [17], which qualitatively compares pollen counts and corona displays. To seek a possible application of advanced photographic analysis on pollen corona displays, we extracted an angular brightness and chromaticity profile from photographic images of corona displays taken by a digital still camera during cryptomeria pollination season, and confirmed that the particles responsible are cryptomeria pollen. We then reproduced the observed pollen corona by Mie theory to retrieve the column number density, assuming the presence of additional aerosols. This article describes the method and results of the attempt. In Section 2 we present the observed pollen corona and the meteorological conditions for the display. Section 3 introduces the method and the results of particle radius estimation, and Section 4 describes the method of reproduction and estimation of pollen column number density. A summary is presented in Section 5. 2. Photographic Image and Meteorological Conditions
A bright solar corona was observed at Tohoku University Aobayama Campus, Sendai, Japan on
Fig. 1. Photographic image of the solar corona taken at Tohoku University on April 12, 2012, 10:05 a.m. The white circles denote the angle from the center of the Sun (1°, 2°, and 3°). The black triangle in the lower half is a building used to shade the camera from direct solar radiation.
April 12, 2012. The corona was first confirmed at 8:40 a.m. (Japan Standard Time: JST, UTC + 9) and continuously displayed until the last observation at 1:00 p.m. The corona disappeared after the last observation as cirrus spread. Figure 1 is a photographic image taken with a digital camera (Fujifilm F70EXR) at 10:05 a.m. The corona at that time was bright enough to be observed with the naked eye. To prevent direct solar radiation straying into the camera, the solar disk was concealed by a building (black triangle in the lower portion of the image). The rings of the solar corona were concentric and circularly symmetric, with the first red ring located at 2° from the Sun. The weather at the time the photograph was taken was clear and no clouds were seen. The absence of clouds through most of the solar corona display is supported by the attenuated backscattering coefficient, which was measured by a lidar system installed and owned by the National Institute for Environmental Studies (NIES) (Fig. 2). The Mie lidar have outputs of 20 mJ (532 nm) and 30 mJ (1064 nm), a repetition rate of 10 Hz, temporal resolution of 5 min, and spatial resolution of 60 m. The lidar measures the attenuated backscattering coefficient at two wavelengths and the depolarization ratio at 532 nm. Figure 2 shows the attenuated backscattering coefficients as measured by lidar; no cloud is seen up to 15 km altitude, with a layer of aerosol up to 1500 m. Figure 3 shows precipitation events, temperature, and relative humidity at the Aobayama Campus, and the pollen count at Tohoku University Hospital (2.4 km from the Aobayama Campus). An increased pollen count precedes the appearance of the corona; the pollen count started increasing at 7:00 a.m., peaking at 500 particles∕m at around noon. A gradual decrease in humidity also preceded the pollen count; humidity started decreasing at 3:00 a.m. on April 12, and reached 20% at 2:00 p.m. The 1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B13
Fig. 2. Attenuated backscatter coefficient observed by the NIES Mie lidar. The arrow shows the period of corona display. High clouds appeared in the last hour of the period of the corona display. The aerosol layer is confined to within 1500 m of the ground.
temperature continuously increased during the morning, providing favorable conditions for pollen release (a temperature increase accompanied by a humidity decrease). Rain events were observed from 2:00 p.m. on April 11 to 2:00 a.m. on April 12. Because of the rain, the atmosphere contained less aerosol than normal at the time of the solar corona display. The optical thickness was 0.1, according to independent sunphotometer observation. Therefore, the particles causing the corona were newly released on April 12 or advected by the northwesterly wind. Figure 4 shows the wind direction at 9:00 a.m. and nearby vegetation. Upwind of the observation point (point a) is a wide patchy distribution of cryptomeria forest, including several dense forests (e.g., area b, 35 km from point a). From these meteorological conditions and the vegetation distribution, it is highly likely that the particles responsible for the corona display were cryptomeria pollen.
Fig. 3. Precipitation, temperature, relative humidity, and pollen count on the day of the pollen corona display and the preceding day. Multiple rain events were recorded from the afternoon of April 11 to the early hours of April 12. The humidity decreased and the pollen count increased before the corona display was confirmed. Pollen count data are from the Pollen Observing System, Ministry of the Environment. B14
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
Fig. 4. Vegetation distribution and surface winds at 9:00 a.m., April 12. The red area is cryptomeria forest, green is other vegetation. Urbanized areas and agricultural land are white. Patchy cryptomeria forest spreads over 100 km upwind of the observation point a, including dense forest in area b (35 km from point a). Vegetation data are based on the Vegetation Survey of the 2nd–5th National Survey of the Natural Environment, Ministry of the Environment (1998). The wind direction was obtained from the Japan Meteorological Agency.
3. Estimating Particle Size A. Atmospheric Model
The two-layer model illustrated in Fig. 5 is utilized to simulate radiance near the forward scattering, based on the finding that the vertical profile of attenuation backscattering coefficient measured by the NIES lidar (Fig. 2) is confined to below 1500 m above ground level. The stratospheric ozone absorption is simulated in the first (upper) model layer, and scattering by aerosols is treated in the second (lower) layer. Reflection from the ground is assumed to be negligible. Scattering by atmospheric molecules is confined in the second layer. The two-layer model employed here is equivalent to a single-slab isolated homogeneouslayer model with stratospheric ozone absorption. Column ozone density is fixed at 300 DU. To calculate the single scattering albedo and extinction cross section by Mie theory [18], all aerosol particles are assumed to be spherical.
Fig. 5. Conceptual diagram of the two-layer model used to estimate the particle size and reproduce the corona display. The model is equivalent to an isolated homogeneous-layer model with ozone absorption.
When the pollen corona image (Fig. 1) was taken, the solar elevation was as high as 54°, and the aerosol optical thickness was 0.1 according to an independent sunphotometer observation. These conditions justify the use of a plane-parallel model with single-scattering approximation in the second layer. The solution of the radiative transfer equation under single-scattering approximation is obtained analytically as the following equations. Below, Iτ; −μ; ϕ is the radiance at the bottom of the two-layer model in the direction of −μ; ϕ, when the top layer is lit by the irradiance, F 0 . The direction of direct radiation is −μ0 ; ϕ0 , the optical thickness of the second layer is τ, and μ cos θ for a zenith angle, θ. The scattering angle, Θ, is determined by a function of μ, ϕ, μ0 , and ϕ0 . ( Iτ; −μ; ϕ
ϖ μ0 F 0 PΘ exp − μτ − exp 4π μ−μ0 ϖ τF 0 PΘ exp − μτ0 4π μ0
− μτ0
μ ≠ μ0 μ μ0
.
(1) The single-scattering albedo, ϖ, and the phase function, PΘ P−μ; ϕ; μ; ϕ, are calculated by using Mie theory. The phase function is normalized to satisfy the condition 1 4π
Z 4π
P−μ; ϕ; μ; ϕdΩ 1:
(2)
In some prior studies, the finite size of the solar disk, ϕs 0.26°, is considered [11,13] because the apparent size of the Sun blurs the vivid color separation. To accomplish this, we used the convoluted ˆ phase function, PΘ, instead of the original phase function, PΘ. They are related by Eq. (3), Z Θϕ s ˆ PΘ0 wΘ0 ; ΘdΘ0 ; (3) PΘ Θ−ϕs
wΘ0 ; Θ
is a weighting function that depends where on the scattering angle, Θ. The weighting function is obtained by differentiating Eq. (4) with respect to ϕ, as follows: Z
l Θ2 l lΘ − 2; 2G wΘ ; ΘdΘ F ϕs ϕ ϕs ϕs Θ−ϕs ϕ
0
0
(4)
where functions Fx, Gx, and l are defined as Fx Gx
l
arcsinx ϕ2 − ϕ2s − Θ2 ≤ 0 π − arcsinx ϕ2 − ϕ2s − Θ2 > 0;
(5)
ϕ2s − ϕ2 − Θ2 ≤ 0 ϕ2s − ϕ2 − Θ2 > 0;
(6)
arcsinx π − arcsinx
1∕2 1 2 2 1 2 ϕs ϕ − ϕ ϕ2s − Θ2 2 : Θ 4
In the derivation of Eqs. (4)–(7), we assumed that the solar disk has uniform brightness. This treatment of the finite size of the Sun more rigorously take into account the effect of two-dimensional convolution than Eq. (10) in [13]. The wavelengths at which radiances are calculated spans from 380 to 780 μm in 5 μm steps, summing up to 61. To simulate the angular distribution and chromaticity (when the whole image is not needed), Eq. (1) is used because the camera’s angle of view is limited and meridional variation is averaged out. B. Estimating Radius by x -Value Profile
To estimate the size of the particles responsible for the corona, the angular distribution of brightness is analyzed. From diffraction theory, the angular sizes of dark rings θ1 and θ2 at wavelength, λ, provide a good approximation of the particle radius, r (μm), by the following equation: 8 j1s < j1s λ π 2 sin θ1 π 1.21967 . (8) r j j2s : 2s λ 2.23313 π 2 sin θ2 π Several different wavelengths, λ, have been used in prior studies, including the older λ 0.570 μm and the widely used λ 0.490 μm. The use of λ 0.570 μm is justified by the assumption that an apparent red, bright ring appears at the green, dark ring [1]. A shorter wavelength of λ 0.490 μm is popularly used [2,13,19] because it can correct the effect of anomalous diffraction [2], which is the difference between the diffraction pattern and the Mie scattering pattern. Shaw and Neiman [3] took a different approach, using a bright ring. The conventional methods mentioned above all relate a specific wavelength to an apparent coronal ring color, leaving the subjectivity to the investigator’s color recognition. We developed an alternative method that objectively evaluates digital image RGB color profiles without relating monochromatic wavelengths to an apparent color. The method converts radiances calculated by the two-layer model at 61 different wavelengths to the x-value in CIE-1931 color space, which is defined by the International Commission on Illumination (CIE). The x-value minimum calculated from the model and the digital image are compared to estimate the particle radius. x-value is a measure of redness, large when the color is close to red and small when the color is close to blue. The x-value is calculated from the radiance, Iλ, at wavelength, λ, by the following equations: Z X D
λmin
Z (7)
YD
λmax
λmax λmin
¯ IλXλdλ;
(9)
¯ IλYλdλ;
(10)
1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B15
Z ZD
x
λmax
λmin
¯ IλZλdλ;
X ; X Y Z
(11)
(12)
where D is a constant that depends on ISO speed, ¯ aperture, and shutter speed. The functions Xλ, ¯ ¯ Yλ, and Zλ are color matching functions of a standard chromatic observer, defined in CIE-1931. The original digital image is initially corrected for vignetting with an empirical correction function, then the RGB values are linearized with respect to radiance. The linear RGB values obtained Rlin ; Glin ; Blin are converted to X, Y, Z values by the matrix operation in Eq. (13) to acquire an observed x-value by Eq. (12), 0
1 1 0 X Rlin @ Y A M @ Glin A; Blin Z
(13)
where the 3 × 3 real matrix M is determined by the response function of the camera and the color matching functions of CIE-1931. We assumed that the camera response function follows the ideal response function over sRGB color space. This assumption can make the result slightly uncertain, as individual cameras may have different response functions to overcome technical difficulties and to satisfy consumer needs. Finally, the x-values calculated for individual pixels are aggregated as an angular profile of x-values. The angular profile of x-values is compared to the angular profile calculated by the two-layer model, varying radii from 10 to 20 μm. The minima of each x-value profile are identified and compared to estimate the particle radius. C.
particles; this dependence is a shared feature with the minima estimated by diffraction theory (line without circles). By utilizing the dependence, particle radius can be estimated from x-minima position. The wavelength assumed for the diffraction theory is λ 0.6 μm. The particle radius is estimated from Fig. 7 to be about 14.5 μm, which is consistent with the range of cryptomeria pollen radii (12–16 μm) as reported in biological studies by Ueno [15]. This consistency supports the judgment that the particles responsible for the corona display are cryptomeria pollen. Table 1 summarizes the radii estimated by Mie theory and diffraction theory. The radii estimated from the first minimum and the second minimum are closer when Mie theory is used. Although more scenes should be used to test the claim, this result suggests the use of Mie theory can possibly
Results
An angular profile of x-values is constructed from the photograph presented in Fig. 1 by the method described in the previous section. The profile is displayed in Fig. 6 along with the Y-value profile, which is a measure of brightness, and the y-value profile, which is a measure of greenness. The brightness sharply decreases with increased scattering angle, supporting the existence of a strong forward-scattering peak. The x-value oscillates with period 1.2°, implying that the separated corona color is periodic. The x-value minima reside at 1.50° and 2.67°, with 95% confidence intervals of 1.48°–1.51° and 2.62°–2.71°, respectively. The confidence intervals are estimated from the noise level of the x-profile. The x-value minima calculated by the two-layer model, with the particle size distribution being a δ-function, are shown in Fig. 7 (line with circles and triangles). The scattering angles at which x-values reach minima depend on the radii of the B16
Fig. 6. x-value profile constructed from the photograph (Fig. 1). The solid line is the Y-value (brightness), the dashed line is the x-value (redness), and the dotted line is the y-value (greenness).
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
Fig. 7. Angle at which x-values reach minima as a function of particle radius. The first minima corresponds to the solid line and the second minima to the dashed line. The lines with circles and triangles are estimated by the two-layer model together with Mie theory, and the lines without markers are estimated by diffraction theory (λ 0.6 μm). The crosshair represents the observed angle with 95% confidence intervals estimated from the noise level.
Table 1.
Theory Mie Diffraction
Particle Radius, as Estimated by Mie Theory and Diffraction Theory
Table 2.
Aerosol
First Minimum
Second Minimum
14.5 0.2 μm 14.0 μm
14.2 0.3 μm 14.4 μm
contribute to more accurate estimation of pollen radius. Note that this comparison became possible by directly processing the digital image without intervention of human color recognition. In the following sections, we assume that the radius of pollen is 14.5 μm because the first minimum is more distinctive than the second. 4. Reproduction of the Corona Display
To further investigate aerosols that induce corona displays, the observed corona is numerically reproduced with the effect of other aerosols and molecular Rayleigh scattering. The two-layer model described in Section 3.A is utilized in the reproduction. A.
Determining the Parameters
To simulate angular profiles of brightness and chromaticity, an aerosol particle size distribution is needed as an input to the two-layer model. Limiting the aerosol to pollen only, the simulated brightness and chromaticity did not well match with observation; we therefore added water soluble aerosol and dust aerosol to minimize the discrepancy. The addition of new aerosol species to the particle size distribution required nine particle size parameters. The range of those parameters is narrowed down by the observed Y-value angular profile. In total, 105 parameter sets were prepared to simulate the corona display. The particle size distribution of the aerosol is assumed to be a trimodal lognormal distribution. An individual lognormal mode is formulated as Eq. (14),
dN N 1 lnr − lnrmod 2 exp − : dln r 2π1∕2 lnσ 2 lnσ (14) The parameters N, σ, and rmod are number density, the geometric standard deviation, and the mode radius, respectively. The geometric standard deviation determines the particle size distribution width, and the mode radius is related to effective radius reff by Eq. (15), reff rmod exp
5 ln σ2 : 2
Preparing Size Distribution Parameters
(15)
The photograph of the corona display captures the scattering angle from 0° to 4°. The brightness profile in this range flattens as the particle size decreases. Thus, the smaller the particle size, the less particle size distribution information is contained. With regard to this size dependence of information content,
Water Soluble
Dust
Pollen
Complex OPAC-WASO OPAC-DUST OPAC-WASO refractive index Effective radius 0.078 μm 2, 3, 4 μm 14.5 μm Geometric standard 2.24 2.6, 1.8, 1.3 1.01, 1.1, 1.2, deviation 1.3, 1.4 Linear Linear Number density 0.05, 0.1, 0.5, inversion inversion 1, 7, 10, 20, 50 × 1012 m−2
different approaches for each aerosol species are taken to prepare the size distribution parameter. For water soluble aerosols, the effective radius and the geometric standard deviation (σ) are specified from OPAC [20]. For dust aerosols, the effective radius and the geometric standard deviation are selected with the aid of the Y-value angular profile (described in Section 4.B). The effective radius of pollen is fixed, but multiple values are prescribed for the geometric standard deviation. Eight densities of water soluble aerosols are prescribed to maintain consistency with the observed Y-value profile, and the densities of dust and pollen are determined by linear inversion, as described in Section 4.C. Table 2 summarizes the approach used to preparing the parameter set. B. Size Distribution of Dust Particles by Y -Value
The phase function for dust aerosols at scattering angles 2°–3° is sensitive to the particle size distribution. Thus, the possible particle size distribution can be limited to several pairs of parameters by using the observed Y-value profile. To identify such pairs, multiple dust size distributions are applied to the two-layer model, and the calculated Y-value profile is compared to the observed profile. Seven effective radii (reff 0.5, 1, 2, 3, 4, 5, 8 μm) and 18 geometric standard deviations (σ 1.1, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6) are prescribed, forming 126 parameter sets. The Y-value profile calculated by the twolayer model is normalized so that Yˆ 1 at 2.18° and Yˆ 0 at 3°. This normalization is intended to remove contributions from pollen and water soluble aerosols, which have nearly constant phase functions. Also, the effect of dust number density is masked out. The observed Y-value is normalized in the same way. The similarity of the modeled Y-value profile to the observed profile is evaluated by the residual square sum (RSS) between the two, as RSS
N 1X ˆ i − Yˆ obs Θi 2 ; YΘ N i1
(16)
ˆ i is the angular distribution of the calcuwhere YΘ ˆ lated Y-value, and Yˆ obs Θi is the angular distribuˆ tion of the Y-value obtained from the photograph. 1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B17
P Θ
Pw ϖ w τ w Pd ϖ d τ d Pp ϖ p τ p PR τ R : (20) τ ϖ
In the equations above, subscript w denotes water soluble aerosol, d denotes dust, p denotes pollen, and R denotes molecular Rayleigh scattering. To simplify Eq. (17), the optical thicknesses of water soluble particles is assumed to be sufficiently larger than that of dust and pollen (τd τp ≪ τw ). The assumption justifies the approximation,
τ exp − μ0 Fig. 8. Variation of Yˆ curve shape with different geometric standard deviation, σ (effective radius 2 μm). The curve is upward convex when the σ is small and downward convex when the σ is large.
The calculated RSS determines the possible combination of an effective radius and the geometric standard deviation (σ). In the scattering angle range ˆ from 2.18° to 3°, the shape of Y-value is determined by σ. Figure 8 shows that, at a fixed effective radius, the curve is upward convex when the σ is small and downward convex when the σ is large. The systematic change of the curve shape implies that the geometric standard deviation σ with smallest RSS best reproduces the observation for a given effective radius. Some effective radii have no RSS minimum in the prescribed σ range, which implies the best σ is out of the probable range for dust aerosols. Such effective radii are removed from the analysis. The best pairs of effective radius and σ obtained by the method described above are reff ;σ 2 μm;2.6, (3 μm, 1.8), (4 μm, 1.3). These three dust aerosol particle size distributions are used in the numerical reproduction of the pollen corona display. C.
Number Density of Dust and Pollen
The number densities of dust and pollen are determined by a simple linear inversion. Assuming three aerosol species—dust, pollen, and water soluble—the observed radiance is a function of abundance and the scattering angle. The solution with single scattering approximation is applicable [Eq. (1)],
F0 τ I ϖ τ P Θ exp − : 4πμ0 μ0
(17)
Here, ϖ , τ , P Θ are total single scattering albedo, total optical thickness, and total phase function, respectively. They are calculated from individual single scattering albedo, φx , optical thicknesses, τx , and phase functions, Px Θ, as follows: τ τw τd τp τR ; ϖ B18
ϖ w τw ϖ d τd ϖ p τp τR ; τ
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
(18)
(19)
τw τR ≈ exp − : μ0
(21)
Substituting the expression into Eq. (17) gives I
F0 τ τR exp − w 4πμ0 μ0 ×Pw Θϖ w τw Pd Θϖ d τd Pp Θϖ p τp PR ΘτR :
(22)
Three new functions, aΘ, bΘ, cΘ, defined as F0 τw τR exp − Pd Θϖ d ; aΘ 4πμ0 μ0
(23)
F0 τ τR exp − w Pp Θϖ p ; 4πμ0 μ0
(24)
bΘ
F0 τw τR Pw Θϖ w τw PR ΘτR ; cΘ exp − 4πμ0 μ0 (25) make it clear that the reduced radiance, I − c, is a linear combination of functions depending on only Θ; namely, aΘ and bΘ, and optical thicknesses (τd , τp ). IΘ − cΘ aΘτd bΘτp :
(26)
Since the Y-value is defined as a linear combination of radiance, I, at different wavelengths [see Eq. (10)], the Y-value is also a linear combination of functions depending only on Θ and optical thicknesses. For a set of discretized Θ (Θ1 ; Θ2 ; ; ΘN ), the Y-value is expressed in a matrix form as 0 B B B B @
1
0 AΘ1 C B C B AΘ2 CB . C @ . . A AΘN YΘN − CΘN YΘ1 − CΘ1 YΘ2 − CΘ2 .. .
1 BΘ1 BΘ2 C C τd .. C τ : . A p
BΘN
(27)
The matrix equation is viewed as simultaneous linear equations with respect to τd and τp . The solutions for τd and τp are easily obtained because functions AΘ, BΘ, CΘ can be computed by the two-layer model with appropriate particle size distribution, and YΘ is already obtained from the photographic analysis. The number density is determined from the solved optical thicknesses τd ; τp via linear approximation. D.
Results
Various corona appearances are reproduced as a result of allowed combinations of particle size distributions. In total, 105 different coronas are simulated from seven size distributions for a water soluble aerosol, three for dust, and five for pollen. The best combination of size distributions is determined by the RSS of the Y-value [see Eq. (16)] from scattering angle 0.88° to 3.06°. The limitation on the range of scattering angle is due to two reasons: the observed brightness saturates at a scattering angle smaller than 0.88°, and vignetting is substantial at scattering angles larger than 3.05°. Figure 9 shows the best reproduction of the observed corona display, and the Y-value profile is presented in Fig. 10. Comparison of the upper and lower halves of Fig. 9 reveals that the chromaticity and radii of the rings are well reproduced. The overall blueness of the photograph is likely to be because the response function of the camera deviated from the ideal sRGB response function as a result of the automatic white point correction. The particle size distribution used for the reproduction shown in Fig. 9 is as follows: the number density of the water soluble aerosol is 1 × 1013 m−2 , that of dust is 1.2 × 108 m−2 , the effective radius of dust is 3 μm with a geometric standard deviation
Fig. 9. Reproduced corona (lower half of the image) and observed corona (the same as Fig. 1, upper half of the image). The radii of the rings and chromaticity profile are well reproduced. The best reproduction is achieved when the number density of the water soluble aerosol is 1 × 1013 particles∕m2 , the effective radius of dust is 3.0 μm with a geometric standard deviation of 1.8, and the geometric standard deviation of pollen is 1.1.
Fig. 10. Angular profiles of the Y-value reproduced and observed. The discrepancy at scattering angles smaller than 1° is because of saturation, and the discrepancy at scattering angles larger than 3° results from uncertainty in the substantial vignetting correction.
of 1.8, and the number density of pollen is 1.3 × 106 m−2 with an effective radius of 14.5 μm and a geometric standard deviation of 1.1. The optical thicknesses of the water soluble aerosol, dust, and pollen are 0.0225 (83%), 0.0030 (11%), and 0.0015 (5.6%), respectively. The total optical thickness of 0.027 is very small compared with the one observed by the sunphotometer (0.1). The difference in total optical thickness requires precise error estimates for the sunphotometer. Despite the uncertainty in the total optical thickness, the estimated pollen optical thickness (0.0015) is far smaller than that of common aerosol species. The analyzed corona display was surprisingly vivid, probably because other aerosols were removed by wet deposition by rain the previous night; without rain, other abundant aerosol species could have obscured the corona display. Therefore, to estimate the pollen column number density from the pollen corona image, contributions from other aerosol species should be correctly treated in the inversion scheme. The column number density of pollen ranged from 1 × 106 m−2 to 1.4 × 106 m−2 in the best twenty sets of reproductions. Assuming that the pollen is distributed uniformly up to 1500 m above ground level, the volumetric number density for pollen ranges from 670 m−3 to 930 m−3 , which is the same order of magnitude with the number density observed by a pollen counter (Fig. 3). The range is also consistent with the pollen count during the boreal pollen corona display [21]. The relatively stable range of estimated pollen number density, irrespective of the various other aerosol particle size distributions implies the effectiveness of an approach that incorporates other aerosols into the model based on brightness profiles at scattering angles exceeding 2°. This stability is advantageous for the number density estimation. 1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B19
Fig. 11. Changes in corona appearance with different geometric standard deviations. The corona is vivid when the geometric standard deviation (σ) is small, and disappears before the σ reaches 1.2. Difference is most evident within a scattering angle between 1° and 2.5° (in the red box).
From the reproduction experiment, it also became clear that the size distribution of pollen is narrow. Figure 11 shows how vividness of the corona varies with different geometric standard deviation. The corona is vivid when the geometric standard deviation is small (i.e., when the size distribution is narrow) but it becomes faint as the geometric standard deviation increases, vanishing at σ > 1.2. Therefore, a visible corona is displayed when the geometric standard deviation is smaller than 1.2. 5. Summary
The vivid corona display photographed with a digital still camera at Sendai, Japan on April 12, 2012 was attributed to suspended pollen of the cryptomeria (or Japanese cedar, or Sugi; Cryptomeria japonica), based on the meteorological conditions, vegetation, and particle size estimates by Lorenz–Mie theory and single scattering radiative transfer calculations. Utilizing direct simulation of brightness and chromaticity, the size of the particles responsible for the corona display was estimated to be 14.5 μm. In our analysis, camera response functions are assumed to be those derived from ideal CIE 1931 chromatic observer and sRGB specification. Although these results are based on a single image of a corona, it seems likely that more consistent size estimation can be achieved by incorporating complete chromaticity and the Mie scattering into the retrieval scheme, than by applying conventional monochromatic difB20
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
fraction theory. The method we introduced is suitable for quantitative observations of pollen corona since no subjective human color recognition is required. The observed corona display was well reproduced by assuming the presence of a water soluble aerosol and dust, in addition to pollen. The reproduction experiment revealed that the geometric standard deviation of pollen particle size distribution must be smaller than 1.2 to cause a visible corona display. In the reproduction experiment, the range of unknown parameters related to particle size distribution are successfully narrowed down by the observed Y-value profile for scattering angles exceeding 2°. Despite the size distribution parameters being narrowed down and varied, the number concentration of pollen was relatively constant among the twenty best reproductions, and consistent with ground observations and the literature. The optical thickness of pollen is small (0.0015), sharing only 5.6% of total optical thickness even under clean atmospheric conditions, such as after a rainfall event. Still, the vivid corona display carries enough information to estimate particle size and number concentration, by incorporating the contributions of other aerosols that are inferred by brightness profiles at scattering angles exceeding 2°. The photographic analysis of coronas, which is conventionally implemented to estimate cloud particle size, is, thus, also applicable to the particle size and column density estimation of pollen particles. References 1. H. C. van de Hulst, Light Scattering by Small Particles (Wiley, 1957). 2. K. Sassen, “Corona-producing cirrus cloud properties derived from polarization lidar and photographic analyses,” Appl. Opt. 30, 3421–3428 (1991). 3. J. Shaw and P. Neiman, “Coronas and iridescence in mountain wave clouds,” Appl. Opt. 42, 476–485 (2003). 4. J. A. Shaw and N. J. Pust, “Icy wave-cloud lunar corona and cirrus iridescence,” Appl. Opt. 50, F6–F11 (2011). 5. G. Molesini and M. Vannoni, “Atypical features of a lunar corona,” J. Opt. A 8, 423–426 (2006). 6. P. Parviainen, C. Bohren, and V. Makela, “Vertical elliptic coronas caused by pollen,” Appl. Opt. 33, 4548–4551 (1994). 7. E. Trankle and B. Mielke, “Simulation and analysis of pollen coronas,” Appl. Opt. 33, 4552–4562 (1994). 8. W. Schneider and M. Vollmer, “Experimental simulations of pollen coronas,” Appl. Opt. 44, 5746–5753 (2005). 9. F. Mims, “Solar corona caused by juniper pollen in Texas,” Appl. Opt. 37, 1486–1488 (1998). 10. P. Laven, “Simulation of rainbows, coronas, and glories by use of Mie theory,” Appl. Opt. 42, 436–444 (2003). 11. S. Gedzelman and J. Lock, “Simulating coronas in color,” Appl. Opt. 42, 497–504 (2003). 12. S. Gedzelman, “Simulating halos and coronas in their atmospheric environment,” Appl. Opt. 47, H157–H166 (2008). 13. J. Lock and L. Yang, “Mie Theory Model of the Corona,” Appl. Opt. 30, 3408–3414 (1991). 14. R. Hayashi, F. Hyodo, J. Urabe, and H. Takahara, “Changes in the accumulation rate of Cryptomeria japonica pollen in sediments from Lake Biwa during the last 100 years,” Palynol. Soc. Jpn. 58, 5–17 (2012) (in Japanese). 15. J. Ueno, Kafun-gaku Kenkyu (Study of Palynology) (KazamaShobo, 1987) (in Japanese).
16. N. Inamura, Y. Shibahara, M. Ishigaki, T. Takasaka, and Y. Sato, “The influence of meteorological factors upon the pollen counts of Cryptomeria japonica,” Nippon Jibiinkoka Gakkai Kaiho 91, 907–914 (1988) (in Japanese). 17. Y. Takahashi, S. Kawashima, and J. Pikki, “Relationship between corona caused by pollen and the numbers of airborne pollen,” Jpn. J. Palynol. 45, 153–157 (1999) (in Japanese). 18. C. Bohren and D. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, 1983).
19. K. Sassen, G. Mace, J. Hallett, and M. Poellot, “Coronaproducing ice clouds: a case study of a cold mid-latitude cirrus layer,” Appl. Opt. 37, 1477–1485 (1998). 20. M. Hess, P. Koepke, and I. Schult, “Optical properties of aerosols and clouds: the software package OPAC,” Bull. Am. Meteorol. Soc. 79, 831–844 (1998). 21. K. Sassen, “Boreal tree pollen sensed by polarization lidar: depolarizing biogenic chaff,” Geophys. Res. Lett. 35, L18810 (2008).
1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B21
Center of Atmospheric and Oceanic Studies, Graduate School of Science, Tohoku University, 6-3 Aramaki-aza-Aoba, Sendai 980-8578, Japan 2
Currently at Department of Atmospheric Sciences, Texas A&M University, MS 3150 TAMU, College Station, Texas 77843-3150, USA *Corresponding author: [email protected] Received 2 July 2014; revised 13 August 2014; accepted 14 August 2014; posted 15 August 2014 (Doc. ID 215247); published 2 October 2014
Brightness and chromaticity profiles were extracted from a vivid solar corona image taken with a digital camera in Sendai, Japan, to compare with a radiative transfer simulation applying Lorenz–Mie theory and single-scattering approximation. The comparison revealed suspended particles having a narrow particle size distribution peaking at radius 14.5 μm. Presumably, pollen of an indigenous coniferous tree, the cryptomeria (Cryptomeria japonica), is responsible for the corona display. The extracted brightness and chromaticity profiles are reproduced well by assuming the presence of a water soluble aerosol and dust in addition to the pollen. We find that photographic analysis of corona displays, similar to that used to measure cloud particle size, is applicable to estimating pollen particle size distribution and column number density. © 2014 Optical Society of America OCIS codes: (010.1290) Atmospheric optics; (010.1110) Aerosols; (010.1310) Atmospheric scattering; (010.1690) Color. http://dx.doi.org/10.1364/AO.54.000B12
1. Introduction
Coronas are iridescent concentric rings around the Sun or the Moon. This atmospheric optical display consists of multiple color-separated rings with apparent radii of a few degrees and is present when particles with radii of a few micrometers are uniformly suspended. As the radii of suspended particles can be estimated from the apparent size of a corona using diffraction theory [1], the estimated particle size is often reported with corona displays due to cloud particles and pollens. Observational approaches to coronas widely utilize photographic analysis in estimating particle radii [2–4]. Among them, Sassen (1991) [2] presents the largest number of samples. Some observed cloud 1559-128X/15/040B12-10$15.00/0 © 2015 Optical Society of America B12
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
coronas have noncircular shapes [3,5], implying effects of cloud heterogeneity or phase transition. Noncircular coronas can be caused by pollen. In 1992, a Finnish amateur astronomers’ network reported vertical elliptical coronas [6]. A simulation of elliptical coronas that uses the scattering theory of spheroids showed that nonspherical pollens can induce noncircular or partly brightened coronas, due to aerodynamic preferences of falling particle orientation [7]. The production of irregular coronal displays from particle orientation was confirmed by experimental simulation [8]. The shape and the size of pollen corona rings depend on the pollen species. Coronas caused by birch, pine, and juniper pollens have been reported [9]. Numerical simulations of pollen corona displays are limited because of difficulties in the calculation of scattering properties of nonspherical particles and radiative transfer with oriented particles.
However, cloud corona displays have become increasingly accessible over the past decade. As a result of increased computing power, scattering calculations for spherical particles by using Lorenz–Mie theory (Mie theory, below) are now possible even on personal computers. Laven showed that fast Mie theory can simulate rainbows, coronas, and glories [10]. Recent simulations of cloud corona displays realistically include background Rayleigh scattering, multiple scattering by clouds, the curvature effect of the Earth, and the effect of ozone absorption [11,12]. This is a major improvement from earlier works such as [13], which calculated scattering patterns by water cloud particles only. Pollen often causes serious allergic symptoms, and is, thus, a public health concern. In Japan, mature coniferous forests of cypress and cryptomeria (or Japanese cedar, or Sugi; Cryptomeria japonica), have been increasing since the 1960s as a result of governmental afforestation policies and the decline of Japan’s forestry industry in later years [14]. These trees cause allergy symptoms among about a third of the population of Japan. Cryptomeria pollen, mainly released from March to April of each year, has been well documented regarding size, allergens, and preferential weather conditions for release. Cryptomeria pollen is nearly spherical, except for a small bump on its top, with a radius ranging from 12 to 16 μm [15]. The pollen count in Sendai, Japan reportedly increases when the temperature increases and humidity concurrently decreases, especially after rain [16]. Despite atmospheric pollen having a great impact on human activities, there have been only limited optical analyses of pollen corona displays. To the best of the authors’ knowledge, the only photographic analysis of pollen coronas in Japan is that of Takahashi et al. [17], which qualitatively compares pollen counts and corona displays. To seek a possible application of advanced photographic analysis on pollen corona displays, we extracted an angular brightness and chromaticity profile from photographic images of corona displays taken by a digital still camera during cryptomeria pollination season, and confirmed that the particles responsible are cryptomeria pollen. We then reproduced the observed pollen corona by Mie theory to retrieve the column number density, assuming the presence of additional aerosols. This article describes the method and results of the attempt. In Section 2 we present the observed pollen corona and the meteorological conditions for the display. Section 3 introduces the method and the results of particle radius estimation, and Section 4 describes the method of reproduction and estimation of pollen column number density. A summary is presented in Section 5. 2. Photographic Image and Meteorological Conditions
A bright solar corona was observed at Tohoku University Aobayama Campus, Sendai, Japan on
Fig. 1. Photographic image of the solar corona taken at Tohoku University on April 12, 2012, 10:05 a.m. The white circles denote the angle from the center of the Sun (1°, 2°, and 3°). The black triangle in the lower half is a building used to shade the camera from direct solar radiation.
April 12, 2012. The corona was first confirmed at 8:40 a.m. (Japan Standard Time: JST, UTC + 9) and continuously displayed until the last observation at 1:00 p.m. The corona disappeared after the last observation as cirrus spread. Figure 1 is a photographic image taken with a digital camera (Fujifilm F70EXR) at 10:05 a.m. The corona at that time was bright enough to be observed with the naked eye. To prevent direct solar radiation straying into the camera, the solar disk was concealed by a building (black triangle in the lower portion of the image). The rings of the solar corona were concentric and circularly symmetric, with the first red ring located at 2° from the Sun. The weather at the time the photograph was taken was clear and no clouds were seen. The absence of clouds through most of the solar corona display is supported by the attenuated backscattering coefficient, which was measured by a lidar system installed and owned by the National Institute for Environmental Studies (NIES) (Fig. 2). The Mie lidar have outputs of 20 mJ (532 nm) and 30 mJ (1064 nm), a repetition rate of 10 Hz, temporal resolution of 5 min, and spatial resolution of 60 m. The lidar measures the attenuated backscattering coefficient at two wavelengths and the depolarization ratio at 532 nm. Figure 2 shows the attenuated backscattering coefficients as measured by lidar; no cloud is seen up to 15 km altitude, with a layer of aerosol up to 1500 m. Figure 3 shows precipitation events, temperature, and relative humidity at the Aobayama Campus, and the pollen count at Tohoku University Hospital (2.4 km from the Aobayama Campus). An increased pollen count precedes the appearance of the corona; the pollen count started increasing at 7:00 a.m., peaking at 500 particles∕m at around noon. A gradual decrease in humidity also preceded the pollen count; humidity started decreasing at 3:00 a.m. on April 12, and reached 20% at 2:00 p.m. The 1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B13
Fig. 2. Attenuated backscatter coefficient observed by the NIES Mie lidar. The arrow shows the period of corona display. High clouds appeared in the last hour of the period of the corona display. The aerosol layer is confined to within 1500 m of the ground.
temperature continuously increased during the morning, providing favorable conditions for pollen release (a temperature increase accompanied by a humidity decrease). Rain events were observed from 2:00 p.m. on April 11 to 2:00 a.m. on April 12. Because of the rain, the atmosphere contained less aerosol than normal at the time of the solar corona display. The optical thickness was 0.1, according to independent sunphotometer observation. Therefore, the particles causing the corona were newly released on April 12 or advected by the northwesterly wind. Figure 4 shows the wind direction at 9:00 a.m. and nearby vegetation. Upwind of the observation point (point a) is a wide patchy distribution of cryptomeria forest, including several dense forests (e.g., area b, 35 km from point a). From these meteorological conditions and the vegetation distribution, it is highly likely that the particles responsible for the corona display were cryptomeria pollen.
Fig. 3. Precipitation, temperature, relative humidity, and pollen count on the day of the pollen corona display and the preceding day. Multiple rain events were recorded from the afternoon of April 11 to the early hours of April 12. The humidity decreased and the pollen count increased before the corona display was confirmed. Pollen count data are from the Pollen Observing System, Ministry of the Environment. B14
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
Fig. 4. Vegetation distribution and surface winds at 9:00 a.m., April 12. The red area is cryptomeria forest, green is other vegetation. Urbanized areas and agricultural land are white. Patchy cryptomeria forest spreads over 100 km upwind of the observation point a, including dense forest in area b (35 km from point a). Vegetation data are based on the Vegetation Survey of the 2nd–5th National Survey of the Natural Environment, Ministry of the Environment (1998). The wind direction was obtained from the Japan Meteorological Agency.
3. Estimating Particle Size A. Atmospheric Model
The two-layer model illustrated in Fig. 5 is utilized to simulate radiance near the forward scattering, based on the finding that the vertical profile of attenuation backscattering coefficient measured by the NIES lidar (Fig. 2) is confined to below 1500 m above ground level. The stratospheric ozone absorption is simulated in the first (upper) model layer, and scattering by aerosols is treated in the second (lower) layer. Reflection from the ground is assumed to be negligible. Scattering by atmospheric molecules is confined in the second layer. The two-layer model employed here is equivalent to a single-slab isolated homogeneouslayer model with stratospheric ozone absorption. Column ozone density is fixed at 300 DU. To calculate the single scattering albedo and extinction cross section by Mie theory [18], all aerosol particles are assumed to be spherical.
Fig. 5. Conceptual diagram of the two-layer model used to estimate the particle size and reproduce the corona display. The model is equivalent to an isolated homogeneous-layer model with ozone absorption.
When the pollen corona image (Fig. 1) was taken, the solar elevation was as high as 54°, and the aerosol optical thickness was 0.1 according to an independent sunphotometer observation. These conditions justify the use of a plane-parallel model with single-scattering approximation in the second layer. The solution of the radiative transfer equation under single-scattering approximation is obtained analytically as the following equations. Below, Iτ; −μ; ϕ is the radiance at the bottom of the two-layer model in the direction of −μ; ϕ, when the top layer is lit by the irradiance, F 0 . The direction of direct radiation is −μ0 ; ϕ0 , the optical thickness of the second layer is τ, and μ cos θ for a zenith angle, θ. The scattering angle, Θ, is determined by a function of μ, ϕ, μ0 , and ϕ0 . ( Iτ; −μ; ϕ
ϖ μ0 F 0 PΘ exp − μτ − exp 4π μ−μ0 ϖ τF 0 PΘ exp − μτ0 4π μ0
− μτ0
μ ≠ μ0 μ μ0
.
(1) The single-scattering albedo, ϖ, and the phase function, PΘ P−μ; ϕ; μ; ϕ, are calculated by using Mie theory. The phase function is normalized to satisfy the condition 1 4π
Z 4π
P−μ; ϕ; μ; ϕdΩ 1:
(2)
In some prior studies, the finite size of the solar disk, ϕs 0.26°, is considered [11,13] because the apparent size of the Sun blurs the vivid color separation. To accomplish this, we used the convoluted ˆ phase function, PΘ, instead of the original phase function, PΘ. They are related by Eq. (3), Z Θϕ s ˆ PΘ0 wΘ0 ; ΘdΘ0 ; (3) PΘ Θ−ϕs
wΘ0 ; Θ
is a weighting function that depends where on the scattering angle, Θ. The weighting function is obtained by differentiating Eq. (4) with respect to ϕ, as follows: Z
l Θ2 l lΘ − 2; 2G wΘ ; ΘdΘ F ϕs ϕ ϕs ϕs Θ−ϕs ϕ
0
0
(4)
where functions Fx, Gx, and l are defined as Fx Gx
l
arcsinx ϕ2 − ϕ2s − Θ2 ≤ 0 π − arcsinx ϕ2 − ϕ2s − Θ2 > 0;
(5)
ϕ2s − ϕ2 − Θ2 ≤ 0 ϕ2s − ϕ2 − Θ2 > 0;
(6)
arcsinx π − arcsinx
1∕2 1 2 2 1 2 ϕs ϕ − ϕ ϕ2s − Θ2 2 : Θ 4
In the derivation of Eqs. (4)–(7), we assumed that the solar disk has uniform brightness. This treatment of the finite size of the Sun more rigorously take into account the effect of two-dimensional convolution than Eq. (10) in [13]. The wavelengths at which radiances are calculated spans from 380 to 780 μm in 5 μm steps, summing up to 61. To simulate the angular distribution and chromaticity (when the whole image is not needed), Eq. (1) is used because the camera’s angle of view is limited and meridional variation is averaged out. B. Estimating Radius by x -Value Profile
To estimate the size of the particles responsible for the corona, the angular distribution of brightness is analyzed. From diffraction theory, the angular sizes of dark rings θ1 and θ2 at wavelength, λ, provide a good approximation of the particle radius, r (μm), by the following equation: 8 j1s < j1s λ π 2 sin θ1 π 1.21967 . (8) r j j2s : 2s λ 2.23313 π 2 sin θ2 π Several different wavelengths, λ, have been used in prior studies, including the older λ 0.570 μm and the widely used λ 0.490 μm. The use of λ 0.570 μm is justified by the assumption that an apparent red, bright ring appears at the green, dark ring [1]. A shorter wavelength of λ 0.490 μm is popularly used [2,13,19] because it can correct the effect of anomalous diffraction [2], which is the difference between the diffraction pattern and the Mie scattering pattern. Shaw and Neiman [3] took a different approach, using a bright ring. The conventional methods mentioned above all relate a specific wavelength to an apparent coronal ring color, leaving the subjectivity to the investigator’s color recognition. We developed an alternative method that objectively evaluates digital image RGB color profiles without relating monochromatic wavelengths to an apparent color. The method converts radiances calculated by the two-layer model at 61 different wavelengths to the x-value in CIE-1931 color space, which is defined by the International Commission on Illumination (CIE). The x-value minimum calculated from the model and the digital image are compared to estimate the particle radius. x-value is a measure of redness, large when the color is close to red and small when the color is close to blue. The x-value is calculated from the radiance, Iλ, at wavelength, λ, by the following equations: Z X D
λmin
Z (7)
YD
λmax
λmax λmin
¯ IλXλdλ;
(9)
¯ IλYλdλ;
(10)
1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B15
Z ZD
x
λmax
λmin
¯ IλZλdλ;
X ; X Y Z
(11)
(12)
where D is a constant that depends on ISO speed, ¯ aperture, and shutter speed. The functions Xλ, ¯ ¯ Yλ, and Zλ are color matching functions of a standard chromatic observer, defined in CIE-1931. The original digital image is initially corrected for vignetting with an empirical correction function, then the RGB values are linearized with respect to radiance. The linear RGB values obtained Rlin ; Glin ; Blin are converted to X, Y, Z values by the matrix operation in Eq. (13) to acquire an observed x-value by Eq. (12), 0
1 1 0 X Rlin @ Y A M @ Glin A; Blin Z
(13)
where the 3 × 3 real matrix M is determined by the response function of the camera and the color matching functions of CIE-1931. We assumed that the camera response function follows the ideal response function over sRGB color space. This assumption can make the result slightly uncertain, as individual cameras may have different response functions to overcome technical difficulties and to satisfy consumer needs. Finally, the x-values calculated for individual pixels are aggregated as an angular profile of x-values. The angular profile of x-values is compared to the angular profile calculated by the two-layer model, varying radii from 10 to 20 μm. The minima of each x-value profile are identified and compared to estimate the particle radius. C.
particles; this dependence is a shared feature with the minima estimated by diffraction theory (line without circles). By utilizing the dependence, particle radius can be estimated from x-minima position. The wavelength assumed for the diffraction theory is λ 0.6 μm. The particle radius is estimated from Fig. 7 to be about 14.5 μm, which is consistent with the range of cryptomeria pollen radii (12–16 μm) as reported in biological studies by Ueno [15]. This consistency supports the judgment that the particles responsible for the corona display are cryptomeria pollen. Table 1 summarizes the radii estimated by Mie theory and diffraction theory. The radii estimated from the first minimum and the second minimum are closer when Mie theory is used. Although more scenes should be used to test the claim, this result suggests the use of Mie theory can possibly
Results
An angular profile of x-values is constructed from the photograph presented in Fig. 1 by the method described in the previous section. The profile is displayed in Fig. 6 along with the Y-value profile, which is a measure of brightness, and the y-value profile, which is a measure of greenness. The brightness sharply decreases with increased scattering angle, supporting the existence of a strong forward-scattering peak. The x-value oscillates with period 1.2°, implying that the separated corona color is periodic. The x-value minima reside at 1.50° and 2.67°, with 95% confidence intervals of 1.48°–1.51° and 2.62°–2.71°, respectively. The confidence intervals are estimated from the noise level of the x-profile. The x-value minima calculated by the two-layer model, with the particle size distribution being a δ-function, are shown in Fig. 7 (line with circles and triangles). The scattering angles at which x-values reach minima depend on the radii of the B16
Fig. 6. x-value profile constructed from the photograph (Fig. 1). The solid line is the Y-value (brightness), the dashed line is the x-value (redness), and the dotted line is the y-value (greenness).
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
Fig. 7. Angle at which x-values reach minima as a function of particle radius. The first minima corresponds to the solid line and the second minima to the dashed line. The lines with circles and triangles are estimated by the two-layer model together with Mie theory, and the lines without markers are estimated by diffraction theory (λ 0.6 μm). The crosshair represents the observed angle with 95% confidence intervals estimated from the noise level.
Table 1.
Theory Mie Diffraction
Particle Radius, as Estimated by Mie Theory and Diffraction Theory
Table 2.
Aerosol
First Minimum
Second Minimum
14.5 0.2 μm 14.0 μm
14.2 0.3 μm 14.4 μm
contribute to more accurate estimation of pollen radius. Note that this comparison became possible by directly processing the digital image without intervention of human color recognition. In the following sections, we assume that the radius of pollen is 14.5 μm because the first minimum is more distinctive than the second. 4. Reproduction of the Corona Display
To further investigate aerosols that induce corona displays, the observed corona is numerically reproduced with the effect of other aerosols and molecular Rayleigh scattering. The two-layer model described in Section 3.A is utilized in the reproduction. A.
Determining the Parameters
To simulate angular profiles of brightness and chromaticity, an aerosol particle size distribution is needed as an input to the two-layer model. Limiting the aerosol to pollen only, the simulated brightness and chromaticity did not well match with observation; we therefore added water soluble aerosol and dust aerosol to minimize the discrepancy. The addition of new aerosol species to the particle size distribution required nine particle size parameters. The range of those parameters is narrowed down by the observed Y-value angular profile. In total, 105 parameter sets were prepared to simulate the corona display. The particle size distribution of the aerosol is assumed to be a trimodal lognormal distribution. An individual lognormal mode is formulated as Eq. (14),
dN N 1 lnr − lnrmod 2 exp − : dln r 2π1∕2 lnσ 2 lnσ (14) The parameters N, σ, and rmod are number density, the geometric standard deviation, and the mode radius, respectively. The geometric standard deviation determines the particle size distribution width, and the mode radius is related to effective radius reff by Eq. (15), reff rmod exp
5 ln σ2 : 2
Preparing Size Distribution Parameters
(15)
The photograph of the corona display captures the scattering angle from 0° to 4°. The brightness profile in this range flattens as the particle size decreases. Thus, the smaller the particle size, the less particle size distribution information is contained. With regard to this size dependence of information content,
Water Soluble
Dust
Pollen
Complex OPAC-WASO OPAC-DUST OPAC-WASO refractive index Effective radius 0.078 μm 2, 3, 4 μm 14.5 μm Geometric standard 2.24 2.6, 1.8, 1.3 1.01, 1.1, 1.2, deviation 1.3, 1.4 Linear Linear Number density 0.05, 0.1, 0.5, inversion inversion 1, 7, 10, 20, 50 × 1012 m−2
different approaches for each aerosol species are taken to prepare the size distribution parameter. For water soluble aerosols, the effective radius and the geometric standard deviation (σ) are specified from OPAC [20]. For dust aerosols, the effective radius and the geometric standard deviation are selected with the aid of the Y-value angular profile (described in Section 4.B). The effective radius of pollen is fixed, but multiple values are prescribed for the geometric standard deviation. Eight densities of water soluble aerosols are prescribed to maintain consistency with the observed Y-value profile, and the densities of dust and pollen are determined by linear inversion, as described in Section 4.C. Table 2 summarizes the approach used to preparing the parameter set. B. Size Distribution of Dust Particles by Y -Value
The phase function for dust aerosols at scattering angles 2°–3° is sensitive to the particle size distribution. Thus, the possible particle size distribution can be limited to several pairs of parameters by using the observed Y-value profile. To identify such pairs, multiple dust size distributions are applied to the two-layer model, and the calculated Y-value profile is compared to the observed profile. Seven effective radii (reff 0.5, 1, 2, 3, 4, 5, 8 μm) and 18 geometric standard deviations (σ 1.1, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6) are prescribed, forming 126 parameter sets. The Y-value profile calculated by the twolayer model is normalized so that Yˆ 1 at 2.18° and Yˆ 0 at 3°. This normalization is intended to remove contributions from pollen and water soluble aerosols, which have nearly constant phase functions. Also, the effect of dust number density is masked out. The observed Y-value is normalized in the same way. The similarity of the modeled Y-value profile to the observed profile is evaluated by the residual square sum (RSS) between the two, as RSS
N 1X ˆ i − Yˆ obs Θi 2 ; YΘ N i1
(16)
ˆ i is the angular distribution of the calcuwhere YΘ ˆ lated Y-value, and Yˆ obs Θi is the angular distribuˆ tion of the Y-value obtained from the photograph. 1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B17
P Θ
Pw ϖ w τ w Pd ϖ d τ d Pp ϖ p τ p PR τ R : (20) τ ϖ
In the equations above, subscript w denotes water soluble aerosol, d denotes dust, p denotes pollen, and R denotes molecular Rayleigh scattering. To simplify Eq. (17), the optical thicknesses of water soluble particles is assumed to be sufficiently larger than that of dust and pollen (τd τp ≪ τw ). The assumption justifies the approximation,
τ exp − μ0 Fig. 8. Variation of Yˆ curve shape with different geometric standard deviation, σ (effective radius 2 μm). The curve is upward convex when the σ is small and downward convex when the σ is large.
The calculated RSS determines the possible combination of an effective radius and the geometric standard deviation (σ). In the scattering angle range ˆ from 2.18° to 3°, the shape of Y-value is determined by σ. Figure 8 shows that, at a fixed effective radius, the curve is upward convex when the σ is small and downward convex when the σ is large. The systematic change of the curve shape implies that the geometric standard deviation σ with smallest RSS best reproduces the observation for a given effective radius. Some effective radii have no RSS minimum in the prescribed σ range, which implies the best σ is out of the probable range for dust aerosols. Such effective radii are removed from the analysis. The best pairs of effective radius and σ obtained by the method described above are reff ;σ 2 μm;2.6, (3 μm, 1.8), (4 μm, 1.3). These three dust aerosol particle size distributions are used in the numerical reproduction of the pollen corona display. C.
Number Density of Dust and Pollen
The number densities of dust and pollen are determined by a simple linear inversion. Assuming three aerosol species—dust, pollen, and water soluble—the observed radiance is a function of abundance and the scattering angle. The solution with single scattering approximation is applicable [Eq. (1)],
F0 τ I ϖ τ P Θ exp − : 4πμ0 μ0
(17)
Here, ϖ , τ , P Θ are total single scattering albedo, total optical thickness, and total phase function, respectively. They are calculated from individual single scattering albedo, φx , optical thicknesses, τx , and phase functions, Px Θ, as follows: τ τw τd τp τR ; ϖ B18
ϖ w τw ϖ d τd ϖ p τp τR ; τ
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
(18)
(19)
τw τR ≈ exp − : μ0
(21)
Substituting the expression into Eq. (17) gives I
F0 τ τR exp − w 4πμ0 μ0 ×Pw Θϖ w τw Pd Θϖ d τd Pp Θϖ p τp PR ΘτR :
(22)
Three new functions, aΘ, bΘ, cΘ, defined as F0 τw τR exp − Pd Θϖ d ; aΘ 4πμ0 μ0
(23)
F0 τ τR exp − w Pp Θϖ p ; 4πμ0 μ0
(24)
bΘ
F0 τw τR Pw Θϖ w τw PR ΘτR ; cΘ exp − 4πμ0 μ0 (25) make it clear that the reduced radiance, I − c, is a linear combination of functions depending on only Θ; namely, aΘ and bΘ, and optical thicknesses (τd , τp ). IΘ − cΘ aΘτd bΘτp :
(26)
Since the Y-value is defined as a linear combination of radiance, I, at different wavelengths [see Eq. (10)], the Y-value is also a linear combination of functions depending only on Θ and optical thicknesses. For a set of discretized Θ (Θ1 ; Θ2 ; ; ΘN ), the Y-value is expressed in a matrix form as 0 B B B B @
1
0 AΘ1 C B C B AΘ2 CB . C @ . . A AΘN YΘN − CΘN YΘ1 − CΘ1 YΘ2 − CΘ2 .. .
1 BΘ1 BΘ2 C C τd .. C τ : . A p
BΘN
(27)
The matrix equation is viewed as simultaneous linear equations with respect to τd and τp . The solutions for τd and τp are easily obtained because functions AΘ, BΘ, CΘ can be computed by the two-layer model with appropriate particle size distribution, and YΘ is already obtained from the photographic analysis. The number density is determined from the solved optical thicknesses τd ; τp via linear approximation. D.
Results
Various corona appearances are reproduced as a result of allowed combinations of particle size distributions. In total, 105 different coronas are simulated from seven size distributions for a water soluble aerosol, three for dust, and five for pollen. The best combination of size distributions is determined by the RSS of the Y-value [see Eq. (16)] from scattering angle 0.88° to 3.06°. The limitation on the range of scattering angle is due to two reasons: the observed brightness saturates at a scattering angle smaller than 0.88°, and vignetting is substantial at scattering angles larger than 3.05°. Figure 9 shows the best reproduction of the observed corona display, and the Y-value profile is presented in Fig. 10. Comparison of the upper and lower halves of Fig. 9 reveals that the chromaticity and radii of the rings are well reproduced. The overall blueness of the photograph is likely to be because the response function of the camera deviated from the ideal sRGB response function as a result of the automatic white point correction. The particle size distribution used for the reproduction shown in Fig. 9 is as follows: the number density of the water soluble aerosol is 1 × 1013 m−2 , that of dust is 1.2 × 108 m−2 , the effective radius of dust is 3 μm with a geometric standard deviation
Fig. 9. Reproduced corona (lower half of the image) and observed corona (the same as Fig. 1, upper half of the image). The radii of the rings and chromaticity profile are well reproduced. The best reproduction is achieved when the number density of the water soluble aerosol is 1 × 1013 particles∕m2 , the effective radius of dust is 3.0 μm with a geometric standard deviation of 1.8, and the geometric standard deviation of pollen is 1.1.
Fig. 10. Angular profiles of the Y-value reproduced and observed. The discrepancy at scattering angles smaller than 1° is because of saturation, and the discrepancy at scattering angles larger than 3° results from uncertainty in the substantial vignetting correction.
of 1.8, and the number density of pollen is 1.3 × 106 m−2 with an effective radius of 14.5 μm and a geometric standard deviation of 1.1. The optical thicknesses of the water soluble aerosol, dust, and pollen are 0.0225 (83%), 0.0030 (11%), and 0.0015 (5.6%), respectively. The total optical thickness of 0.027 is very small compared with the one observed by the sunphotometer (0.1). The difference in total optical thickness requires precise error estimates for the sunphotometer. Despite the uncertainty in the total optical thickness, the estimated pollen optical thickness (0.0015) is far smaller than that of common aerosol species. The analyzed corona display was surprisingly vivid, probably because other aerosols were removed by wet deposition by rain the previous night; without rain, other abundant aerosol species could have obscured the corona display. Therefore, to estimate the pollen column number density from the pollen corona image, contributions from other aerosol species should be correctly treated in the inversion scheme. The column number density of pollen ranged from 1 × 106 m−2 to 1.4 × 106 m−2 in the best twenty sets of reproductions. Assuming that the pollen is distributed uniformly up to 1500 m above ground level, the volumetric number density for pollen ranges from 670 m−3 to 930 m−3 , which is the same order of magnitude with the number density observed by a pollen counter (Fig. 3). The range is also consistent with the pollen count during the boreal pollen corona display [21]. The relatively stable range of estimated pollen number density, irrespective of the various other aerosol particle size distributions implies the effectiveness of an approach that incorporates other aerosols into the model based on brightness profiles at scattering angles exceeding 2°. This stability is advantageous for the number density estimation. 1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B19
Fig. 11. Changes in corona appearance with different geometric standard deviations. The corona is vivid when the geometric standard deviation (σ) is small, and disappears before the σ reaches 1.2. Difference is most evident within a scattering angle between 1° and 2.5° (in the red box).
From the reproduction experiment, it also became clear that the size distribution of pollen is narrow. Figure 11 shows how vividness of the corona varies with different geometric standard deviation. The corona is vivid when the geometric standard deviation is small (i.e., when the size distribution is narrow) but it becomes faint as the geometric standard deviation increases, vanishing at σ > 1.2. Therefore, a visible corona is displayed when the geometric standard deviation is smaller than 1.2. 5. Summary
The vivid corona display photographed with a digital still camera at Sendai, Japan on April 12, 2012 was attributed to suspended pollen of the cryptomeria (or Japanese cedar, or Sugi; Cryptomeria japonica), based on the meteorological conditions, vegetation, and particle size estimates by Lorenz–Mie theory and single scattering radiative transfer calculations. Utilizing direct simulation of brightness and chromaticity, the size of the particles responsible for the corona display was estimated to be 14.5 μm. In our analysis, camera response functions are assumed to be those derived from ideal CIE 1931 chromatic observer and sRGB specification. Although these results are based on a single image of a corona, it seems likely that more consistent size estimation can be achieved by incorporating complete chromaticity and the Mie scattering into the retrieval scheme, than by applying conventional monochromatic difB20
APPLIED OPTICS / Vol. 54, No. 4 / 1 February 2015
fraction theory. The method we introduced is suitable for quantitative observations of pollen corona since no subjective human color recognition is required. The observed corona display was well reproduced by assuming the presence of a water soluble aerosol and dust, in addition to pollen. The reproduction experiment revealed that the geometric standard deviation of pollen particle size distribution must be smaller than 1.2 to cause a visible corona display. In the reproduction experiment, the range of unknown parameters related to particle size distribution are successfully narrowed down by the observed Y-value profile for scattering angles exceeding 2°. Despite the size distribution parameters being narrowed down and varied, the number concentration of pollen was relatively constant among the twenty best reproductions, and consistent with ground observations and the literature. The optical thickness of pollen is small (0.0015), sharing only 5.6% of total optical thickness even under clean atmospheric conditions, such as after a rainfall event. Still, the vivid corona display carries enough information to estimate particle size and number concentration, by incorporating the contributions of other aerosols that are inferred by brightness profiles at scattering angles exceeding 2°. The photographic analysis of coronas, which is conventionally implemented to estimate cloud particle size, is, thus, also applicable to the particle size and column density estimation of pollen particles. References 1. H. C. van de Hulst, Light Scattering by Small Particles (Wiley, 1957). 2. K. Sassen, “Corona-producing cirrus cloud properties derived from polarization lidar and photographic analyses,” Appl. Opt. 30, 3421–3428 (1991). 3. J. Shaw and P. Neiman, “Coronas and iridescence in mountain wave clouds,” Appl. Opt. 42, 476–485 (2003). 4. J. A. Shaw and N. J. Pust, “Icy wave-cloud lunar corona and cirrus iridescence,” Appl. Opt. 50, F6–F11 (2011). 5. G. Molesini and M. Vannoni, “Atypical features of a lunar corona,” J. Opt. A 8, 423–426 (2006). 6. P. Parviainen, C. Bohren, and V. Makela, “Vertical elliptic coronas caused by pollen,” Appl. Opt. 33, 4548–4551 (1994). 7. E. Trankle and B. Mielke, “Simulation and analysis of pollen coronas,” Appl. Opt. 33, 4552–4562 (1994). 8. W. Schneider and M. Vollmer, “Experimental simulations of pollen coronas,” Appl. Opt. 44, 5746–5753 (2005). 9. F. Mims, “Solar corona caused by juniper pollen in Texas,” Appl. Opt. 37, 1486–1488 (1998). 10. P. Laven, “Simulation of rainbows, coronas, and glories by use of Mie theory,” Appl. Opt. 42, 436–444 (2003). 11. S. Gedzelman and J. Lock, “Simulating coronas in color,” Appl. Opt. 42, 497–504 (2003). 12. S. Gedzelman, “Simulating halos and coronas in their atmospheric environment,” Appl. Opt. 47, H157–H166 (2008). 13. J. Lock and L. Yang, “Mie Theory Model of the Corona,” Appl. Opt. 30, 3408–3414 (1991). 14. R. Hayashi, F. Hyodo, J. Urabe, and H. Takahara, “Changes in the accumulation rate of Cryptomeria japonica pollen in sediments from Lake Biwa during the last 100 years,” Palynol. Soc. Jpn. 58, 5–17 (2012) (in Japanese). 15. J. Ueno, Kafun-gaku Kenkyu (Study of Palynology) (KazamaShobo, 1987) (in Japanese).
16. N. Inamura, Y. Shibahara, M. Ishigaki, T. Takasaka, and Y. Sato, “The influence of meteorological factors upon the pollen counts of Cryptomeria japonica,” Nippon Jibiinkoka Gakkai Kaiho 91, 907–914 (1988) (in Japanese). 17. Y. Takahashi, S. Kawashima, and J. Pikki, “Relationship between corona caused by pollen and the numbers of airborne pollen,” Jpn. J. Palynol. 45, 153–157 (1999) (in Japanese). 18. C. Bohren and D. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, 1983).
19. K. Sassen, G. Mace, J. Hallett, and M. Poellot, “Coronaproducing ice clouds: a case study of a cold mid-latitude cirrus layer,” Appl. Opt. 37, 1477–1485 (1998). 20. M. Hess, P. Koepke, and I. Schult, “Optical properties of aerosols and clouds: the software package OPAC,” Bull. Am. Meteorol. Soc. 79, 831–844 (1998). 21. K. Sassen, “Boreal tree pollen sensed by polarization lidar: depolarizing biogenic chaff,” Geophys. Res. Lett. 35, L18810 (2008).
1 February 2015 / Vol. 54, No. 4 / APPLIED OPTICS
B21
