Materials and methods

Virtual particles with various physical characteristics of bioaerosol particles were simulated based on equation (19) of Miki (Reference Miki2020) with the addition of the vertical wind velocity, W, and the neglect of the gravitational effect of the sun. In the simulation, the motion of a particle was assumed to be governed by gravity, buoyancy, friction from the air and photophoretic force, which as discussed by Wainwright et al. (Reference Wainwright, Alharbi and Wickramasinghe2006). The photophoretic force results from the temperature difference in the gas next to the particles that results from the difference in solar heating. For micron sized particles exposed to sunlight the photophoretic force can be 20% of the weight of the particle (e.g. Tehranian et al., Reference Tehranian, Giovane and Blum2001). The terminal velocity of a particle (V) is given by:

(1)$$V = \displaystyle{{( \pi D_{{\rm ae}}^3 /6) ( {\rho_{\rm p}-\rho_{\rm f}} ) g-F_{\Delta \alpha } + F_{\Delta {\rm T}}} \over {( ( C_{\rm D}Re_{\rm p}) /24) ( {3\pi \eta D_{{\rm ae}}} ) }}-W, \;$$

where D _ae is the aerodynamic diameter of the particle; W is the vertical wind speed in the stratosphere; r is the radius of the particle; ρ _p is the density of the particle, assumed to be the same as water (1000 kg m⁻³) based on the aerodynamic diameter; ρ _f is the density of air; F _Δα and F _ΔT are the photophoretic forces, which are driven by the difference in the thermal accommodation coefficient inside a particle and the thermal gradient of the particle, respectively; C _D is the drag coefficient; Re _p is the Reynolds number of the particle; and η is the viscosity of air.

The viscosity of air (η) was derived from the Sutherland's formula as:

(2)$$\eta = \left({\displaystyle{T \over {T_0}}} \right)^{3/2} \times \displaystyle{{T_0 + 110.4} \over {T + 110.4}} \times 1.7932 \times 10^{{-}5}, \;$$

where T ₀ is the temperature at sea level and T is the temperature at each altitude, which was calculated as equation (3) using an approximate formula used in the US standard atmosphere (United States Committee on Extension to the Standard Atmosphere, 1976) as first introduced in Miki (Reference Miki2020).

(3)$$T = 30.46 \times \sin \left({\displaystyle{{z-34619.46} \over {28253.7}}\pi } \right) + 238.06.$$

Gryazin and Beresnev (Reference Gryazin and Beresnev2011) reported that in the stratosphere, vertical wind blows at a monthly average speed of approximately 5.0 mm s⁻¹ in upward and downward directions and an annual average speed of approximately 1.0 mm s⁻¹ in upward and downward directions. In this study, wind speeds were set as 5.0 mm s⁻¹ (5.0 × 10⁻³ m s⁻¹) and 1.0 mm s⁻¹ (±1.0 × 10⁻³ m s⁻¹) in upward and downward directions, and 0 mm s⁻¹ because the simulation duration significantly varied depending on the variable settings such as particle size or density.

Reynold's number (Re _p) was derived using the following equation:

(4)$$Re_{\rm p} = \displaystyle{{\rho _{\rm f}( {V + W} ) D_{{\rm ae}}} \over \eta }.$$

The density of air (ρ _f) varied depending on the altitude of the stratosphere, and for each altitude, ρ _f was calculated using equations (5) and (6):

(5)$$\rho _{\rm f} = \displaystyle{M \over {RT}}P, \;$$

(6)$$P( h ) = P_0\exp \left({-\displaystyle{h \over H}} \right), \;$$

where M is the molar mass of the air, R is the gas constant, P is the air pressure at each altitude and H is the scale height of the Earth's atmosphere (8432 m). The approximate formula for the drag coefficient (C _D) changed depending on the amplitude of the Reynold's number of the particle (Re _p) as shown in equations (7) and (8). These were first introduced by Oseen (Reference Oseen1910).

When Re _p  <  1:

(7)$$C_{\rm D} = \displaystyle{{24} \over {Re_{\rm p}}}, \;{\rm \;}$$

and when Re _p  ≥  1:

(8)$$C_{\rm D} = \displaystyle{{24} \over {Re_{\rm p}}}\left({1 + \displaystyle{3 \over {16}}Re_{\rm p}} \right).$$

The two types of photophoretic forces, F _ΔT and F _Δα, were calculated using equations (9)–(14), which have been previously introduced (Rohatschek, Reference Rohatschek1996; Wurm and Krauss, Reference Wurm and Krauss2008). It was assumed that F _ΔT works downwards and F _Δαworks upward (Keith, Reference Keith2010).

(9)$$F_{\Delta {\rm T}} = F_{\rm \ast }\displaystyle{2 \over {( P/p_{\rm \ast }) + ( p_{\rm \ast }/P) }}$$

(10)$$F_{\Delta \alpha } = \displaystyle{1 \over {12c_{\rm m}}}{\rm \;}\displaystyle{{\Delta \alpha } \over \alpha }\pi r^2S\displaystyle{1 \over {{\left({1 + \displaystyle{P \over {\,p_{\rm \ast }}}} \right)}^2}}$$

(11)$$p_{\rm \ast } = D\sqrt {\displaystyle{2 \over \alpha }} \displaystyle{{3T} \over {\pi r}}$$

(12)$$F_{\rm \ast } = D\sqrt {\displaystyle{\alpha \over 2}} \displaystyle{{r^2JS} \over {C_{\rm T}}}$$

(13)$$D = \displaystyle{\pi \over 2}\sqrt {\displaystyle{\pi \over 3}\kappa \displaystyle{{c_{\rm m}\eta } \over T}} $$

(14)$$c_{\rm m} = \sqrt {\displaystyle{8 \over \pi }\displaystyle{{RT} \over M}} , \;$$

where, $F_\ast$ is the maximum force under pressure $p_{\rm \ast }$, J is the asymmetric parameter, C _T is the thermal conductivity, κ is the thermal creep parameter, c _m is the average speed of the atmospheric molecule, α is the thermal accommodation coefficient and Δα is the difference in the thermal accommodation coefficient inside a particle.

When S is the solar radiation intensity at each altitude,

(15)$$S( L ) = S_0( L ) \exp \left({-\displaystyle{{kM} \over R}P_0\mathop \smallint \limits_z^\infty \displaystyle{1 \over T}\exp \left({-\displaystyle{h \over H}} \right){\rm d}h} \right), \;$$

where S ₀ is the solar radiation at wavelength (L) and k is a coefficient.

From equation (15), the relationship between the solar radiation above the atmosphere and that on the Earth's surface is described as follows:

(16)$$\displaystyle{{\Delta S} \over {S_0}} = 1-\exp \left({-35.74 \times \displaystyle{{kM} \over R}P_0} \right), \;$$

where ΔS is the difference between S ₀ and the solar radiation above the atmosphere S _earth.

From equation (16), k is described as follows:

(17)$$k( L ) = \displaystyle{R \over {-35.74 \times MP_0}}\ln \left({\displaystyle{{S_{{\rm earth}}} \over {S_0}}} \right).$$

When E(h) is given as equation (18), the effective irradiance is derived as equation (19).

(18)$$E( h ) = P_0\mathop \smallint \limits_z^\infty \displaystyle{1 \over T}\exp \left({-\displaystyle{h \over H}} \right){\rm d}h$$

(19)$$I_{\rm E} = \displaystyle{1 \over \pi }\mathop \smallint \limits_0^\infty S_0( L ) \exp \left({-\displaystyle{{kM} \over R}E} \right){\rm d}L$$

Equation (19) can be approximated using the piecewise quadrature as equation (20)

(20)$$I_{\rm E} = \displaystyle{1 \over \pi }\mathop \smallint \limits_{L_i=0}^\infty S_0( L_i ) \exp \left({-\displaystyle{{kM} \over R}E} \right){\Delta}L$$

As the actual solar radiation data, ASTM G-173-03, the distributions of power obtained by the National Renewable Energy Laboratory (NREL) were used (https://www.nrel.gov/grid/solar-resource/spectra-am1.5.html 2022/10/27; Fig. 1a). As the radiation data are obtained discretely every 0.5 or 5.0 nm, the effective irradiance was approximately calculated using equation (21):

(21)$$I_{\rm E}\approx \displaystyle{{\left({\mathop \sum \nolimits_{L_i = L_1}^{L_2} S_0( {L_i} ) \exp \left({-\displaystyle{{kM} \over R}E} \right)\Delta L_i + \mathop \sum \nolimits_{L_i = L_1{\rm^{\prime}}}^{L_2{\rm^{\prime}}} S_0( {L_i} ) \exp \left({-\displaystyle{{kM} \over R}E} \right)\Delta L_i} \right)} \over {2\pi }}, \;$$

where L ₁ = 280.0, L ₂ = 3990.0, $L_1^{\prime} = 280.5$ and $L_2^{\prime} = 3995.0$. The solar irradiance at each altitude was derived as Fig. 1(b).

Fig. 1. (a) Spectrum of solar irradiance at above the atmosphere (S ₀), on the Earth's surface (S _earth), and the difference between S ₀ and S _earth (ΔS). (b) Relationship between the altitude and the total irradiance.

The asymmetric parameter (J) and thermal creep parameter (κ) were set at the values J = 0.5 and κ = 1.14 under the assumption that particle surfaces perfectly absorb light, as described in a previous study (Wurm and Krauss, Reference Wurm and Krauss2008). The thermal conductivity (C _T) and thermal accommodation (α and Δα) depend on the characteristics of the particles and ambient air; these parameters of bioaerosol particles have not yet been studied. We assigned 0.3 and 0.9 values to α to account for low and high thermal accommodations, respectively, and 0 and 0.2 to Δα. The values of thermal conductivity (C _T) were set as 0.1 and 10, respectively, to cover a broad range of high and low thermal conductivities. Lastly, in order to analyse bioaerosol particle motion in the stratosphere, the diameters of some representative bioaerosol particles, such as pollen, spores and microorganisms, were selected as aerodynamic diameters in the simulation. For example, an aerodynamic diameter of 10 μm can be applied for small pollen grains because of various primary pollen taxa, such as small Salix pollen. Additionally, the maximum size of possible microorganism found in the stratosphere was found to be 10 μm (Harris et al., Reference Harris, Wickramasinghe and Lloyd2002; Wainwright et al., Reference Wainwright, Alharbi and Wickramasinghe2006), but was 2 μm in another study (Bryan et al., Reference Bryan, Christner and Guzik2019). Spores were represented by a particle with an aerodynamic diameter of 1–3 μm based on a previous morphology study (Reponen et al., Reference Reponen, Grishpun and Conwell2001). Thus, 1, 2, 5 and 10 μm were substituted as the representative aerodynamic diameters in this study.

In summary, every combination of the parameters given below were simulated.

$$\eqalign{& W\colon [ {-}5, \;{-}1, \;0, \;1, \;5] , \;D_{{\rm ae}}\colon [ 1, \;2, \;5, \;10] , \;\cr & \alpha \colon [ 0.3, \;0.9] , \;\Delta \alpha \colon [ 0, \;0.2] , \;C_{\rm T}\colon [ 0.1, \;10] .} $$

The initial altitude was set to 30 000 or 60 000 m, and the terminal velocity was adjusted every 3600 s based on the equations explained above at each altitude.