One-dimensional Microphysics Model of Venusian Clouds from 40 to 100 km: Impact of the Middle-atmosphere Eddy Transport and SOIR Temperature Profile on the Cloud Structure

The model equations used in this study are based on a 1D cloud microphysics model developed by Imamura & Hashimoto (2001). We note, however, that the purpose of using a 1D model is different from Imamura & Hashimoto (2001), who reproduced the equatorial atmosphere by incorporating the upward branch of the Hadley circulation into the 1D model. Here we focus on the globally averaged vertical structure by representing all vertical transport processes, including large-scale advection such as the Hadley circulation, by eddy diffusion, and thus the advection by vertical winds is absent in the equations. The model simulates dynamical and cloud microphysical processes in a 1D column of CO₂ air from 40 to 100 km. As a major update from the previous work by Imamura & Hashimoto (2001), we extended the model top altitude from 70 to 100 km to include the upper haze layer. CN and cloud droplets are modeled with discrete size bins. The size bins are divided into 23 bins between 0.17 and ∼30 μm, in which the volume doubles from one bin to the next. The smallest size bin represents Mode 1 particles, which are assumed to be insoluble CN. Size bins larger than the smallest bin represent droplets, whose liquid composition consists of a binary solution of H₂SO₄ and H₂O. The cloud size distribution is a function of the particle mass m, the altitude z, and the time t. It is calculated based on the 1D continuity equation defined without the advection term (modified after Imamura & Hashimoto 2001):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial C(m,z,t)}{\partial t}=\displaystyle \frac{\partial }{\partial z}\left({K}_{{zz}}(z)\ \rho (z)\ \displaystyle \frac{\partial }{\partial z}\left[\displaystyle \frac{C(m,z,t)}{\rho (z)}\right]\right)\\ -\displaystyle \frac{\partial }{\partial z}\left[{w}_{\mathrm{sed}}(m,z)\ C(m,z,t)\right]-\displaystyle \frac{\partial }{\partial m}\left[G(m,z,t)\ C(m,z,t)\right]\\ +\displaystyle \frac{1}{2}{\int }_{{m}_{\mathrm{cn}}}^{m}{K}_{\mathrm{coag}}\left({m}^{{\prime} },m-{m}^{{\prime} }\right)\ C\left({m}^{{\prime} },z,t\right)\ C\left(m-{m}^{{\prime} },z,t\right)\,{dm}^{\prime} \\ -C(m,z,t){\int }_{{m}_{\mathrm{cn}}}^{{m}_{\max }}{K}_{\mathrm{coag}}\left({m}^{{\prime} },m\right)\ C\left({m}^{{\prime} },z,t\right){dm}^{\prime} ,\end{array}\end{eqnarray} \tag{ 1 }$$

where C(m, z, t) is the cloud particle number density, ρ(z) is the atmospheric density, K_zz (z) is the vertical eddy diffusion coefficient of air, w_sed(m, z) is the sedimentation velocity of cloud particles, G(m, z, t) is the condensational growth rate, ${K}_{\mathrm{coag}}(m,{m}^{{\prime} })$ is the coagulation kernel between cloud particles with mass m and $m^{\prime}$, and m_cn and m_max are the masses of condensational nucleus and particle with maximum size, respectively. The first term on the right-hand side is vertical transport by eddy diffusion. The second term accounts for vertical transport due to particle sedimentation. The third term represents the growth through condensation, and the fourth and fifth terms represent particle size evolution resulting from coagulation.

Brownian diffusion is accounted for in the coagulation process. The CN is treated the same way as in the earlier work (Imamura & Hashimoto 2001). Particles in the smallest size bin act as the CNs. The composition and microphysical properties of CNs are unknown. Therefore, for simplicity, we assume that they do not coagulate with each other. Thermal equilibrium is assumed when calculating the H₂SO₄ mass fraction in the cloud droplets, considering that equilibrium adjustment is achieved by the exchange of H₂O molecules between the droplet surface and the surrounding environment.

Regarding the sedimentation process, while the previous work by Imamura & Hashimoto (2001) used the Stokes velocity, we implemented the Cunningham slip-flow correction as follows (Jacobson 2005):

$$\begin{eqnarray}&&\begin{array}{l}{w}_{\mathrm{sed}}(m,z)=\displaystyle \frac{2r{\left(m\right)}^{2}{\rho }_{p}g}{9\eta }\\ \,\,\times \,\left(1+\mathrm{Kn}(z)\left[A+B\exp \left(-\displaystyle \frac{C}{\mathrm{Kn}(z)}\right)\right]\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where r(m) is the particle radius corresponding to the mass m; ρ_p is the density of the particle; g = 8.7 m s⁻² is the gravitational acceleration; η = 1.5 × 10⁻⁵ kg m⁻¹ s⁻¹ is the viscosity of CO₂ gas; A, B, and C are experimental parameters set to be 1.249, 0.42, and 0.87, respectively, following Kasten (1962); and Kn(z) is the Knudsen number for air. The inclusion of the Cunningham slip-flow correction is crucial for accurately estimating sedimentation velocity in thin air owing to strong dependency of Kn on altitude. This effect, which was not included in earlier work by Imamura & Hashimoto (2001), has limited impact on particle size distribution between 40 and 70 km.

We also consider the chemical production of H₂SO₄ vapor and loss of H₂O vapor through the following net reaction, as included in previous works (Imamura & Hashimoto 2001; McGouldrick 2017; McGouldrick & Barth 2023):

$$\begin{eqnarray}&&\displaystyle \frac{3}{2}{\mathrm{SO}}_{2}+{{\rm{H}}}_{2}{\rm{O}}\to \displaystyle \frac{1}{2}{\rm{S}}+{{\rm{H}}}_{2}{\mathrm{SO}}_{4}.\end{eqnarray} \tag{ 3 }$$

It demonstrates that the production of one H₂SO₄ molecule results in the loss of one H₂O molecule. We updated the production rate of H₂SO₄ using the photochemical model results from Krasnopolsky (2012). The updated production rate is a globally averaged value, consistent with the other settings in the present study. The column integrated production rate is about 5.7 × 10¹⁵ m⁻² s⁻¹, which is approximately half of the rate used in the study by Imamura & Hashimoto (2001). In addition, the peak altitude of the production rate is located at ∼66 km following Krasnopolsky (2012), while it was situated at 61 km in the previous work. We also adopted the CN production rate introduced by Imamura & Hashimoto (2001), who assumed that the CNs are made of elemental sulfur produced by reaction (3) and have the same radius as particles in Mode 1. The production rate of CN, P_cn, is expressed as

$$\begin{eqnarray}&&{P}_{\mathrm{cn}}(z)=\displaystyle \frac{1}{2}\ {P}_{1}(z){\left(\displaystyle \frac{{\rho }_{\mathrm{cn}}}{{M}_{s}}\displaystyle \frac{4}{3}\pi {r}_{\mathrm{cn}}^{3}\right)}^{-1},\end{eqnarray} \tag{ 4 }$$

where P₁(z) is the H₂SO₄ vapor production rate, r_cn and ${\rho }_{\mathrm{cn}}$ are respectively the radius and mass density of CN, and M_s is the molecular mass of elemental sulfur. The photochemical production rate profiles of H₂SO₄ vapor and the CN are shown in Figure 1. It should be noted that if the production rate of CNs matches that of reaction (3), it implies the unavailability of elemental sulfur for the chemistry of this cycle. This is a simplistic assumption and is not corroborated by current findings on the net sulfur production from photochemical models. These models also fail to account for the potential loss of sulfur into droplets, where sulfur primarily reacts with O₂. However, studying the interaction between the sulfur chemical cycle and CN production goes beyond the scope of this paper.

Zoom In Zoom Out Reset image size

Selecting only one CN size bin may overlook the impact of CNs on microphysical processes. Knollenberg & Hunten (1980) found a wide size distribution of CNs that includes significant numbers of particles both larger and smaller than 0.17 μm. Larger CNs might represent a nonnegligible sink of H₂SO₄ and H₂O, while smaller CNs might be less relevant owing to the Kelvin effect, which hinders condensation on smaller particles. Moreover, it is improbable that CNs with radii as large as 0.17 μm originate directly from photochemical processes because the initial CN products are considerably smaller and contain far fewer molecules. These particles grow to sizes around 0.17 μm through nucleation and coagulation over time, leading to a broad size distribution. Thus, the simplified treatment of CNs in the present study may overestimate the number of photochemically produced CNs as large as 0.17 μm. These effects should be investigated by incorporating a CN size distribution in a future study.

As mentioned earlier, the detailed settings for cloud microphysical processes and underlying chemistry are also adapted from Imamura & Hashimoto (2001); we refer the reader to this paper for more details. In the following subsections, we outline the parameters used in case studies and describe the calculation settings.

The nominal eddy diffusion profile is established based on the past estimations (see Figure 2). Eddy diffusion is not an observable quantity. It is rather a conceptual parameterization of atmospheric dynamics over a variety of scales.

Zoom In Zoom Out Reset image size

In detailing the historical perspective, radio scintillation measurements have played an important role. They estimated the contribution of small-scale turbulence to the diffusion coefficient, approximating it to be ∼0.2 m² s⁻¹ at 45 km (Woo et al. 1982) and ∼4 m² s⁻¹ at 60 km (Woo & Ishimaru 1981). To align our estimated eddy diffusion profile with these observations, we interpolate the eddy diffusion coefficient from 0.1 m² s⁻¹ at the bottom boundary (40 km) to 4 m² s⁻¹ at 60 km. Consequently, the interpolated eddy diffusion is ∼0.2 m² s⁻¹ at 45 km, consistent with the aforementioned studies.

For a more nuanced representation, we considered the influence of convective mixing in the cloud due to infrared heating of the cloud base. This phenomenon, evidenced by the neutral static stability observed by entry probes (e.g., Seiff et al. 1980) and the results from PV, VEx, and Akatsuki radio occultation (Kliore & Patel 1982; Tellmann et al. 2009; Imamura et al. 2017), necessitated an adjustment in our model. Following the work by Imamura & Hashimoto (2001), we incorporated an enhanced eddy diffusion coefficient of 250 m² s⁻¹ centered at 53 km.

It is important to note the theoretical insights from Baker et al. (1998, 1999, 2000a, 2000b), which indicated that downwelling plumes within the convective layer may extend into the subcloud region, potentially inducing substantial vertical mixing. This convective mixing process can alter the size distribution of particles, as reported in the studies by Imamura & Hashimoto (2001) and McGouldrick & Toon (2007). While this effect is not extensively examined in the current study, its possible impacts on cloud microphysics will be addressed in Section 3.1.

Between 60 and 70 km, the eddy diffusion profile is held constant at the 60 km value for the convenience of parameter studies. The eddy diffusion above 85 km is determined based on past estimations by von Zahn et al. (1980). The eddy diffusion coefficient was derived by comparing observed gas species profiles with those calculated numerically using a 1D diffusion equation. The formulated coefficient was ∼1.4 × 10⁹[M]^−0.5 m² s⁻¹, yielding a value of ∼10 m² s⁻¹ at 90 km. Thus, we use this value as a typical eddy diffusion coefficient above 85 km. For the layer between 70 and 85 km, the eddy diffusion coefficient is deduced by interpolating the values obtained for the adjacent layers, providing a consistent and gradual transition in the profile.

To investigate the effects of eddy transport processes on particles and gas species, we performed a case study that consisted of four pairs of different eddy diffusion coefficients at 60–70 km and 85–100 km (Cases 1, 2, and 3, in addition to the Nominal Case; see Figure 2). The estimation of eddy diffusion using radio scintillation is based on the assumption that all observed signal fluctuations result from small-scale refractive index structures caused by turbulence; however, the refractive index fluctuation can be due to gravity waves rather than turbulence (Leroy & Ingersoll 1996). Consequently, the eddy diffusion coefficient derived by Woo & Ishimaru (1981) is likely an upper limit. On the other hand, larger-scale eddies, including the mean meridional circulation, can also contribute to vertical transport, and their contributions should be included in the eddy diffusion coefficient in the 1D model. To explore this variability and its impact on particle and gas-phase species distributions, we reduced the eddy diffusion coefficient by a factor of 2 between 60 and 70 km for Cases 2 and 3. Recently, Mahieux et al. (2021) estimated the eddy diffusion coefficient in the mesosphere using a 1D photochemical model and CO and CO₂ profiles obtained by the SOIR instrument on board VEx and He and CO₂ observations by PV. They derived a coefficient of ∼360 m² s⁻¹ from 80 to 100 km, an order of magnitude larger than the value obtained by von Zahn et al. (1980). Although Mahieux et al. (2021) derived a complete data set of the eddy diffusion coefficient ranging from 80 to 130 km, the bottom boundary in their 1D photochemical model was 85 km. Therefore, we use the value of 360 m² s⁻¹ above 85 km for the parameter study (Cases 1 and 3). These four eddy diffusion cases are determined as four patterns composed of pairs of eddy diffusion coefficients for the 60–70 km and 80–100 km altitude ranges, with interpolation between the coefficients for these two ranges.

The background temperature and pressure profiles are taken from the Venus International Reference Atmosphere (VIRA; Seiff et al. 1985) for the Nominal Case and eddy diffusion case studies. We chose the midlatitude profiles to represent the global average of the background atmosphere, as shown in Figure 3. The temperature and pressure profiles are held constant in the calculation.

Zoom In Zoom Out Reset image size

In addition, we defined case studies using a different temperature profile obtained by the VEx observations. Bertaux et al. (2007) found a warm layer in the Venus mesosphere using the SOIR data set, indicating that the temperature above 85 km could be up to 50 K higher than previously reported. Subsequent studies also confirmed the existence of this warm layer during the course of VEx observations (Mahieux et al. 2012, 2015, 2023b). Mahieux et al. (2023b) proposed that the observed significant temperature differences between the VIRA and SOIR data can be attributed to the variations in temperature distribution with respect to local time, since SOIR observations were all confined to the terminator region. The uncertainties on the individual SOIR temperature profiles vary between 20 and 40 K, while the variability is around 50 K at all altitudes. The temperature difference between the VIRA and SOIR measurements could affect particle and gas distributions by changing the saturation vapor pressure. Therefore, we considered the midlatitude temperature profile from Mahieux et al. (2015) to test the sensitivity toward temperature (Figure 3) in Cases 4 and 5. The SOIR temperature profiles are limited to altitudes above 80 km. We integrated these profiles with those from VIRA to establish a comprehensive temperature distribution for all the altitudes of interest in our model. This integration was achieved by seamlessly connecting the SOIR data with the distributions from VIRA. The connection point was carefully chosen at the altitude where the temperature distributions from VIRA and SOIR naturally intersect, ensuring a continuous temperature profile across the entire altitude range under study. We use the same eddy diffusion coefficient for the Nominal Case as for Case 4 as a control experiment with the SOIR temperature profile. On the other hand, Case 5 uses the same eddy diffusion coefficient as in Case 3, as it well replicates the observed particle and gas distributions. The names of each case (Nominal, Case 1, Case 2, Case 3, Case 4, and Case 5), along with their corresponding parameters, are listed in Table 2. It is important to note that the SOIR temperature profiles exhibit large interprofile variability. This variation can impact condensation and evaporation processes, potentially resulting in changes to aerosol size and acidity. However, the study of the microphysical response to such variation is beyond the scope of our current research.

Table 3 shows the initial and boundary conditions in our simulations, which are determined based on past observations. The H₂O vapor VMR is reported to be ∼30 ppm below 40 km (Marcq et al. 2008; Bézard et al. 2011; Chamberlain et al. 2013; Arney et al. 2014; Fedorova et al. 2015). We take this value as the bottom boundary condition for the H₂O vapor. The bottom boundary for the H₂SO₄ vapor VMR is set at 4 ppm, based on previous observations from both Magellan (Kolodner & Steffes 1998) and VEx (Oschlisniok et al. 2021). The in situ observations reported by Knollenberg & Hunten (1980) suggested that the Mode 1 particles extended below the cloud layer. These particles are assumed to be nonvolatile since H₂SO₄ vapor cannot exist in the liquid phase within the subcloud region. Therefore, the particle number density for the smallest radius is set at 4 × 10⁷ m⁻³, based on the number density of the Mode 1 particles from in situ observations. The mixing ratio gradients of cloud particles, H₂SO₄ vapor, and H₂O vapor are set to be zero at the top boundary. The initial conditions for H₂O, H₂SO₄ VMR, and the smallest radius bin are set to match the bottom boundary values across all vertical levels. The aerosol and gas species profiles reach a steady state within ∼50 Venusian years (11,250 Earth days) starting from these initial conditions. This timescale is much longer than the timescales of the physical processes incorporated within the model, such as eddy diffusion, sedimentation, condensation, and coagulation. We performed the simulation for 100 Venusian years (22,500 Earth days) and present the steady-state results in the next section.

Table 3. Boundary Conditions

Altitude	CN Number Density	Droplet Number Density	H2SO4 VMR	H2O VMR
40 km (bottom)	4 × 107 m−3	0 m−3	4 ppm	30 ppm
100 km (top)	$\displaystyle \frac{{{df}}_{\mathrm{CN}}}{{dz}}=0$	$\displaystyle \frac{{{df}}_{\mathrm{droplet}}}{{dz}}=0$	$\displaystyle \frac{{{df}}_{H2{SO}4}}{{dz}}=0$	$\displaystyle \frac{{{df}}_{{\rm{H}}2{\rm{O}}}}{{dz}}=0$

Note. The mixing ratio is expressed as f.

Download table as:  ASCIITypeset image

Figures 4 and 5 show the simulated size distribution in the Nominal Case and Case 1. The particle size distribution can be divided into three regions by its formation processes, being consistent with the previous microphysics studies (e.g., Gao et al. 2014; McGouldrick & Barth 2023): the lower and middle cloud (47–57 km), the upper cloud (57–70 km), and the upper haze (70–100 km) regions. The upper cloud particles are formed through the condensation of photochemically produced H₂SO₄ vapor onto CN, which are also produced photochemically through reaction (3) within the same altitude range. The calculated mode radius for the upper cloud layer is ∼0.7 μm, which corresponds to Mode 2 particles but is somewhat underestimated. The lower and middle clouds exhibit a bimodal size distribution, with peaks at the smallest radius bin (0.17 μm) and ∼3 μm corresponding to Modes 1 and 3, respectively (Figures 5(a) and (c)). Between 51 and 53 km, the Kelvin effect inhibits the condensational growth of submicron particles, allowing larger particles to grow. The evaporation of the submicron particles and the upward diffusion transport of CN from the bottom boundary contribute to the formation of the peak at the smallest bin of 0.17 μm. Conversely, the larger peak forms as a result of condensational growth of the photochemically produced particles descending from the upper cloud region. The condensation process is driven by H₂SO₄-rich air ascending from lower altitudes in the convective layer, where strong eddy diffusion is assumed (see Figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Above 53 km, the H₂SO₄ supersaturation becomes sufficiently high to suppress the Kelvin effect, allowing all particles to undergo condensational growth (Figure 5(b)). Consequently, the smaller peak clearly seen below 53 km altitude eventually merges into the size distribution of the upper cloud. Below 51 km altitude, all particles start to evaporate because of the increasing H₂SO₄ saturation vapor pressure. The evaporation gradually decreases the droplet number density, leading to the near-complete disappearance of droplets around the cloud base at around 47 km altitude.

Figure 5 demonstrates that the size distribution simulated in the present study closely agrees with earlier work by Imamura & Hashimoto (2001). They showed that the size distribution in the lower and middle cloud layers is bimodal at steady state and that a trimodal distribution can form owing to convective mixing. Furthermore, Baker et al. (1998) suggested that the microphysical properties of the clouds could be altered by condensation processes initiated by cold downwelling plumes from neutrally stable regions in the cloud layers. This implies that the absence of an explicit connection with convective activity in the current study may explain why a trimodal distribution is not observed in these cloud layers.

The upper haze features a broad, unimodal size distribution with a mode radius of ∼0.4 μm (Figures 5(e) and (f)). This radius possibly represents Mode 2 particles as observed by Wilquet et al. (2009) or the effective radius calculated by Luginin et al. (2016), assuming a unimodal size distribution. The particle population in the upper haze is maintained by the upward transport of particles from the upper cloud layer via eddy diffusion.

Changes in the eddy diffusion coefficient above 85 km lead to an enhanced upper haze population. It is evident that the upper haze extends vertically more than in the Nominal Case because of efficient eddy diffusion transport (Figure 4(b)). The size distribution closely resembles that of the Nominal Case at 70 km, but the differences between the two cases become increasingly evident at higher altitudes (Figures 5(d), (e), and (f)). The total particle number density of Case 1 is orders of magnitude larger than in the Nominal Case at 90 km, and larger particles are more abundant in Case 1. This underscores the importance of eddy diffusion transport in shaping the upper haze particle population and its microphysical properties.

Figure 6 shows the simulated cloud structure as a function of altitude. To obtain the mass loading, we calculated the density of the H₂SO₄–H₂O solution droplets by interpolating the density of pure H₂O liquid and pure H₂SO₄ liquid based on weight fraction of H₂SO₄ (see Figure 9(a)). For the determination of the extinction coefficient, the open-source software MiePython (Prahl 2023) is used to calculate the Mie scattering extinction efficiency of the aerosols. For the input of Mie scattering calculation, we used the refractive index from Palmer & Williams (1975) as a function of the concentration of the H₂SO₄–H₂O solution. The simulated mass loading in the Nominal Case is basically consistent with the PV in situ observations by Knollenberg & Hunten (1980) (Figure 6(a)). The mass loading is highest in the middle cloud layer owing to rapid cloud formation induced by convective activity, aligning with in situ observations. The simulated value is about a factor of 2 smaller than the observed one in the middle cloud layer. The mass loading in the upper cloud is consistent with the observations, although slightly underestimated by a factor of 2. The in situ observations were conducted at the equator, capturing the photochemical production rates of H₂SO₄ vapor and CN in low latitudes. In contrast, our model uses rates from midlatitudes, which are lower. This likely accounts for the discrepancy between the simulated and observed values.

Zoom In Zoom Out Reset image size

On the other hand, the Nominal Case fails to replicate the observed extinction coefficient profile (Figure 6(b)). In the upper haze layer, the extinction coefficient decreases exponentially with altitude, but the simulated gradient is steeper than observed, resulting in values about two orders of magnitude smaller. The use of increased eddy diffusion coefficient in Cases 1 and 3 reduces the discrepancy with SOIR observations and results in a more extended upper haze layer. It is clearly shown that the cases with efficient eddy diffusion are in good agreement with the observed vertical profiles by Wilquet et al. (2012). This suggests that the recently estimated eddy diffusion coefficient (Mahieux et al. 2021) is suitable for the reproduction of the observed upper haze optical properties.

Figures 7 and 8 show the vertical profiles of condensational gas species. The H₂SO₄ VMR is fixed at 4 ppm at the bottom of the model and takes the maximum value of ∼6 ppm around the cloud bottom at 47 km for all cases (Figure 7(a)). This peak arises from droplet sedimentation in the lower cloud and subsequent evaporation. Between the cloud bottom and 60 km, the H₂SO₄ VMR aligns with the saturation vapor profile. This pattern agrees with observations from the Magellan mission (Kolodner & Steffes 1998) and recent VEx data (Oschlisniok et al. 2021), suggesting H₂SO₄ VMR levels of 5–7 ppm in the subcloud region at midlatitudes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Above 75 km, the simulated H₂SO₄ VMR presents two distinct profiles. The first includes the Nominal Case and Case 2 scenarios with a maximum VMR of ∼5 × 10⁻⁴ ppm; the second includes Cases 1 and 3 with a maximum VMR of ∼5 × 10⁻⁶ ppm in the upper haze layer. The latter cases feature a higher eddy diffusion coefficient above 75 km, leading to a more abundant upper haze particle population. This facilitates condensation onto these particles, as the H₂SO₄ vapor is supersaturated. The peak values for both profile types are lower than the upper limit reported by Sandor et al. (2012), indicating that H₂SO₄ vapor is not a dominant source of sulfur species in the upper haze layer.

The H₂SO₄ VMR is supersaturated above 60 km and reaches its second-highest value at 67 km, due to the photochemical production of H₂SO₄ vapor across all cases. This finding is consistent with Dai et al. (2022), who used a cloud condensation model without detailed microphysics but assumed the same H₂SO₄ production profile as ours. These results indicate that the supersaturation above 60 km is a robust feature across different models when using the photochemical production of H₂SO₄ vapor from the 1D photochemical model of Krasnopolsky (2012). They further showed that if a Gaussian distribution is assumed for H₂SO₄ production, the H₂SO₄ VMR closely follows the saturation VMR in the upper haze layer. This is because the Gaussian distribution exhibits an extremely low value outside the half-width of the peak production altitude, and loss by H₂SO₄ condensation dominates over the photochemical production. The earlier work by Imamura & Hashimoto (2001), which also assumed a Gaussian profile, did not find supersaturation except near the peak production altitude. This sensitivity of the H₂SO₄ VMR profile to the assumed photochemical production profile suggests that it is a critical factor in modeling.

The H₂O vapor VMR in the Nominal Case monotonically decreases from 40 to 67 km and stabilizes at ∼9 ppm above 67 km (Figure 8). This trend aligns qualitatively with VEx observations (Fedorova et al. 2008; Cottini et al. 2012; Fedorova et al. 2016), although the simulated values are slightly higher. The VMR gradient below the cloud base is steeper than in earlier work by Imamura & Hashimoto (2001) and recent microphysics model results reported by McGouldrick (2017) and McGouldrick & Barth (2023). The reason is that we use a lower eddy diffusion coefficient for those altitudes, based on the previous study by Woo et al. (1982). This steeper H₂O VMR gradient helps balance the H₂O vapor supply by upward eddy diffusion flux against loss from chemical reactions. In the lower and middle cloud layers, the gradient is less steep owing to stronger eddy diffusion, which acts to flatten the H₂O vapor VMR profile. The mass fraction of liquid H₂O in the droplets increases (corresponding to a decrease in liquid H₂SO₄ mass fraction, as shown in Figure 9(a)), indicating that the uptake of H₂O vapor by the droplets also contributes to the decreasing H₂O VMR from lower to upper cloud layers. Around 67 km, the H₂O VMR experiences a sharp drop due to efficient chemical loss by H₂SO₄ formation (Equation (3) and Figure 1).

Zoom In Zoom Out Reset image size

The H₂O VMR is highly sensitive to the assumed eddy diffusion profile used in the model (Figure 8). In Cases 2 and 3, the H₂O VMR gradient is steeper between 60 and 70 km than in the Nominal Case and Case 1. Consequently, the H₂O VMR decreases sharply in the upper cloud layer, aligning well with VEx observations (Fedorova et al. 2008; Cottini et al. 2012; Fedorova et al. 2016; Chamberlain et al. 2020) in these two cases. The defining feature of Cases 2 and 3 is the lower eddy diffusion coefficient between 60 and 70 km, suggesting that the mesospheric H₂O vapor levels are influenced by the efficiency of eddy transport at these altitudes, as expected since the chemical loss of H₂O is fixed.

Figure 9(a) shows liquid H₂SO₄ mass fraction in the cloud particles as a function of altitude. Below 50 km, the H₂SO₄ mass fraction exceeds 0.95 and decreases up to 60 km. It then stabilizes at ∼0.75 from 60 to 75 km in the Nominal Case and Case 1 scenarios, while increasing to ∼0.8 in Cases 2 and 3. The simulated values in the cloud layer are consistent with previous observations (Hansen & Hovenier 1974; Barstow et al. 2012; Cottini et al. 2012) in all cases. In particular, Cases 2 and 3 better match the data from Cottini et al. (2012). This vertical profile is also consistent with reported modeling studies (e.g., Krasnopolsky 2015; Stolzenbach et al. 2023), although our results are overestimated compared to Stolzenbach et al. (2023) below 60 km. Stolzenbach et al. (2023) argued that such a difference could be caused by the lower temperature bias of their model. Above 75 km, the H₂SO₄ mass fraction decreases again, reaching a minimum around 90 km before rising again up to the model top altitude for all cases. The decrease of the H₂SO₄ mass fraction between 75 and 90 km is attributed to a sharp temperature decrease, as temperature is the key factor determining thermal equilibrium.

The equilibrium condition is also influenced by the background H₂O VMR, as the H₂SO₄ mass fraction adjusts to equilibrium through the uptake or release of H₂O molecules at the particle surface. Therefore, a lower ambient concentration of H₂O molecules results in reduced uptake by the droplets. As shown in Figure 8, the H₂O VMR above 60 km is about an order of magnitude lower in Cases 2 and 3 compared to the Nominal Case and Case 1. This leads to less H₂O being absorbed into the droplets and, conversely, a higher H₂SO₄ mass fraction within them.

Figure 9(b) shows the H₂O VMR condensed in droplets as a function of altitude. Upon careful examination, the influence of eddy diffusion on the condensed H₂O VMR is noticeable between 60 and 70 km. Specifically, the condensed H₂O VMR is more concentrated at the cloud formation altitude around 66 km in Cases 2 and 3, due to the reduced eddy diffusion coefficient (see Figure 2). Above 70 km, Cases 1 and 3 exhibit higher levels of liquid H₂O compared to the Nominal Case and Case 2. This is attributable to the extended aerosol layer in these two cases facilitated by the efficient eddy transport, leading to the upward transport of liquid H₂O in its liquid phase. The condensed H₂O VMR aligns well with the condensation model of Dai et al. (2022) in Cases 2 and 3. In their study, the liquid H₂O VMR was ∼2 ppm in the upper cloud layer, ∼0.1 ppm at 50 km, and ∼0.2 ppm at 80 km. Note that the eddy diffusion profile they assumed was ∼0.3 m² s⁻¹ at 40 km, ∼1 m² s⁻¹ at 40–50 km, and 1–6 m² s⁻¹ at 60–70 km (see Table 4), resulting in a similar eddy transport timescale to our model.

A twofold difference in the eddy diffusion coefficient in the upper cloud layer leads to an order-of-magnitude difference in H₂O VMR in the mesosphere, as shown in Figure 8. To investigate the sensitivity of H₂O vertical profile with eddy diffusion, we compared the H₂O flux profiles for the Nominal Case and Case 2 in Figure 10; positive fluxes indicate upward transport, and negative fluxes indicate downward transport. Below 60 km altitude, both profiles are identical since the eddy diffusion coefficient is nearly the same at these altitudes. It is noteworthy that the flux balances for the two cases are also similar between 60 and 70 km, even though the eddy diffusion coefficient for Case 2 is twice as large as that of the Nominal Case. Indeed, in a steady-state condition, the upward H₂O vapor flux essentially balances with the sum of the liquid diffusion flux, liquid sedimentation flux, and cumulative chemical loss. However, as demonstrated in Figure 10(c), the difference in liquid flux between the two cases is relatively small, accounting for less than 20% difference in the total sum of the liquid flux and chemical loss. Given that the cumulative chemical loss is a constant parameter in all cases, the sums of the liquid flux and the cumulative chemical loss end up being similar. Consequently, the upward vapor fluxes are also similar, as they are balanced with the downward fluxes.

Zoom In Zoom Out Reset image size

The two different H₂O VMR trends shown in Figure 8 result from the fact that all cases exhibit similar upward H₂O vapor fluxes between 60 and 70 km. In the present study, the species in the gas phase are transported only by eddy diffusion. The H₂O vapor flux due to eddy diffusion ${{\rm{\Phi }}}_{{{\rm{H}}}_{2}{\rm{O}}}$ is expressed as ${{\rm{\Phi }}}_{{{\rm{H}}}_{2}{\rm{O}}}=-{K}_{{zz}}\ {n}_{a}\tfrac{\partial {f}_{{{\rm{H}}}_{2}{\rm{O}}}}{\partial z}\$, where K_zz is the eddy diffusion coefficient, n_a is the number density of air, and ${f}_{{{\rm{H}}}_{2}{\rm{O}}}$ is the H₂O VMR. As shown in Figure 10, ${{\rm{\Phi }}}_{{{\rm{H}}}_{2}{\rm{O}}}$ does not differ significantly in each case and can be regarded as a constant for simplicity. If K_zz is increased by a factor of 2, the H₂O VMR gradient $\tfrac{\partial {f}_{{{\rm{H}}}_{2}{\rm{O}}}}{\partial z}$ must be divided by 2 to yield the same vapor flux ${{\rm{\Phi }}}_{{{\rm{H}}}_{2}{\rm{O}}}$. This is why a twofold difference in the eddy diffusion coefficient leads to a significant change in the extent of H₂O vapor depletion in the mesosphere. Consequently, the steep gradient reduces the H₂O VMR in the upper cloud layer, leading to a H₂O VMR that is about an order of magnitude smaller above 70 km in Cases 2 and 3.

In the present study, the chemical production/loss rate of condensational gas species is fixed. However, as shown in Figure 10, the cumulative H₂O chemical loss also plays a significant role in determining the compensating upward H₂O vapor flux. Namely, if the chemical loss rate is increased and K_zz is fixed, $\tfrac{\partial {f}_{{{\rm{H}}}_{2}{\rm{O}}}}{\partial z}$ should increase to balance with the chemical loss. Given the amplitude of the cumulative chemical loss compared to the liquid flux components, it is expected that a factor of variation of the loss rate may lead to considerable change in $\tfrac{\partial {f}_{{{\rm{H}}}_{2}{\rm{O}}}}{\partial z}$.

The recent cloud condensation model from Dai et al. (2022) showed a similar relationship between eddy diffusion and the H₂O gradient in the upper cloud layer. With the same H₂O chemical loss profile as in our study, the Nominal Case from Dai et al. (2022) used an eddy diffusion coefficient of 1 m² s⁻¹ at 60 km and ∼3 m² s⁻¹ at 65 km. This suggests that the efficiency of eddy transport in their study aligns closely with ours. Their resulting H₂O VMR gradient between 60 and 70 km was steep enough to reduce the H₂O VMR to a few parts per million at 70 km. In contrast, when they used the eddy diffusion profile from Bierson & Zhang (2020), which had a value of ∼4 m² s⁻¹ between 60 and 70 km, they observed a less steep H₂O gradient. The resulting H₂O VMR exceeded 10 ppm at 70 km. Our results also overestimated the H₂O VMR at 70 km when the eddy diffusion is set to 4 m² s⁻¹ between 60 and 70 km. The similarities between our findings and those from Dai et al. (2022) highlight the critical role of the eddy diffusion coefficient in the upper cloud layer. Furthermore, both studies concur that an eddy diffusion of ∼2 m² s⁻¹ serves as a threshold beyond which hygroscopic effects become significant.

The influence of eddy diffusion on the chemical production/loss rates should also be explored. It has been established that the vertical profiles of sulfur species are significantly influenced by eddy diffusion. Krasnopolsky (2012) demonstrated that variations in eddy diffusion coefficients above 55 km can lead to a significant difference in SO₂ concentrations in the upper cloud region. Additionally, Bierson & Zhang (2020) found that the altitude and rate of H₂SO₄ production vary considerably with the eddy diffusion coefficient in the cloud layers. Given the vital role of the chemical production profile in the formation of the upper cloud and the depletion of H₂O, the H₂SO₄ and H₂O fluxes could be altered by the response of sulfur chemistry to eddy diffusion. However, here we focus on the impact of changing eddy diffusion on cloud microphysics and its effect on the H₂SO₄ and H₂O fluxes, and the coupling between photochemistry and cloud microphysics will be addressed in future research.

We examined the sensitivity of aerosol and gas species distributions to variations in the upper temperature profile by comparing the VIRA temperature cases (Nominal Case and Case 3) and the SOIR temperature cases (Cases 4 and 5). The mass loading and extinction coefficient in both scenarios are very similar up to 80 km and in good agreement with the LCPS observations, as shown in Figure 11. Above 80 km, the upper haze is more extended in Cases 3 and 5 because both cases use the updated eddy diffusion by Mahieux et al. (2021), which is consistent with the results shown in Section 3.1. However, we observed that the population of the upper haze particles is slightly reduced when using the temperature profile from SOIR. As will be discussed later, this reduction is due to the aerosol evaporation, which is triggered by the warm layer evident in the SOIR temperature profile.

Zoom In Zoom Out Reset image size

The H₂SO₄ VMR in Case 5, which uses the SOIR temperature profile, closely resembles that of Case 3 up to ∼85 km altitude but diverges significantly from Case 3 above 85 km (see Figure 12(a)). On the other hand, the H₂SO₄ VMR in Case 4 is larger than the Nominal Case, but the amplitude is much less pronounced. The saturation H₂SO₄ VMR calculated using the SOIR profile (Case 5) is orders of magnitude higher than the one calculated using the VIRA profile, as the maximum temperature of the SOIR profile is up to ∼50 K higher than that in the VIRA profile. Due to this elevated saturation VMR, the H₂SO₄ vapor is undersaturated in the upper haze layer above ∼90 km in Cases 4 and 5 (Figure 12(b)). This undersaturation leads to the evaporation of the upper haze particles, in turn increasing the H₂SO₄ VMR in the SOIR temperature cases. However, the deviation from VIRA temperature case is more evident in Case 5 than in Case 4 because the extended upper haze layer supplies H₂SO₄. As a result of the combined effect of enhanced eddy diffusion and higher temperature, the H₂SO₄ VMR aligns with the upper limit suggested by Sandor et al. (2012) in Case 5, suggesting that the mesospheric H₂SO₄ vapor could be detectable under high-temperature conditions.

Zoom In Zoom Out Reset image size

Previous photochemical studies suggested that ∼1 ppm of H₂SO₄ vapor is required to act as a precursor for SO₂ and SO in the upper haze layer (Zhang et al. 2010, 2012). However, the simulated values are orders of magnitude lower than this suggested value, even when using the SOIR temperature (see Figure 12(a)). This discrepancy raises the possibility that sources other than H₂SO₄ aerosol may be needed to provide sufficient sulfur, such as elemental sulfur suggested by Zhang et al. (2012).

The H₂O VMR in Case 4 does not deviate from that obtained in the Nominal Case, as also seen in H₂SO₄ VMR. This is due to the lack of aerosols in the upper haze layer in both cases. However, when Cases 3 and 5 are compared, the H₂O VMR simulated with the SOIR temperature profile diverges notably from that simulated using the VIRA profile, especially above ∼85 km (Figure 13). In Case 5, the H₂O VMR reaches a minimum value at ∼85 km and then increases with altitude. This pattern aligns with the H₂O VMR profile observed by Fedorova et al. (2008) and Chamberlain et al. (2020). The increase in H₂O VMR above 85 km can be attributed to the strong dependence of equilibrium conditions on the background temperature. It is important to note that the diverging altitude can vary depending on how the temperature profile is constructed in the model and may also change when taking into account the interprofile variability discussed in Section 2.3.

Zoom In Zoom Out Reset image size

The equilibrium condition manifests in the liquid H₂SO₄ mass fraction within the droplets. In Cases 4 and 5, there is a significant increase in the H₂SO₄ mass fraction above 85 km, corresponding to the sharp temperature increase in the SOIR profile (Figure 14(a)). This suggests that the droplets cannot maintain a high H₂O content above this altitude. The reduced liquid H₂O VMR in the droplets is corroborated in Case 5 (Figure 14(b)). Specifically, the liquid H₂O VMR experiences a sharp decline above 85 km. This implies that there are more H₂O molecules in gas phase in the SOIR temperature scenario compared to the VIRA temperature case because of the warm layer that leads to H₂O evaporation. Thus, the increasing H₂O VMR in Figure 13 can be attributed to the release of the H₂O molecules from the upper haze particles that have been transported upward from below 85 km.

Zoom In Zoom Out Reset image size

In Section 3.1, we identified a strong dependence of both aerosol and gas species distributions on the eddy diffusion profiles. The aerosol distribution is highly sensitive to the eddy diffusion coefficient above 85 km (Figure 6). Results for Cases 1 and 3, using the eddy diffusion recently estimated by Mahieux et al. (2021), closely match the upper haze distribution observed by SOIR (Wilquet et al. 2012). Furthermore, Cases 2 and 3 using lower eddy diffusion coefficients between 60 and 70 km are consistent with the observed H₂O VMR distribution reported by VEx observations (Fedorova et al. 2008; Cottini et al. 2015; Fedorova et al. 2016; Chamberlain et al. 2020). Given that Case 3 successfully reproduces both the observed gas species and aerosol distributions, the present study suggests that the eddy diffusion should be ∼360 m² s⁻¹ in the upper haze layer and ∼2 m² s⁻¹ in the upper cloud layer, when the chemical production/loss rates of condensational species are fixed.

In our study, as mentioned earlier, the eddy diffusion includes global transport processes such as Hadley-type circulation since our model represents the whole atmosphere with a 1D column. To compare the eddy diffusion coefficient obtained in the present study with those derived from the general circulation model (GCM), we employ timescale estimations. The eddy transport timescale is given by τ = H²/K_zz , where H is the scale height of the atmosphere. On the other hand, the timescale of the global-scale circulation can be estimated by τ = H/w, where w is the typical vertical wind speed. Assuming that eddy transport is mainly driven by the global-scale circulation, the relationship K_zz = H w can be derived. Takagi et al. (2018) obtained a zonally and temporally averaged vertical wind speed magnitude on the order of ∼0.001 m s⁻¹ at 70 km using Atmospheric GCM for the Earth Simulator for Venus (AFES-Venus). With a corresponding scale height of 5 km, these values yield K_zz = 5 m² s⁻¹, which is within the same order of magnitude as the eddy diffusion coefficient assumed between 60 and 70 km in the present study. Navarro et al. (2021) found a zonally and temporally averaged vertical wind speed magnitude on the order of ∼0.1 m s⁻¹ at 90 km using the Institut Pierre-Simon Laplace (IPSL) GCM, and the scale height at the corresponding altitude is ∼3.7 km. These two values give K_zz = 370 m s⁻¹, which is equivalent to the eddy diffusion coefficient assumed above 85 km in our study. A GCM study by Sugimoto et al. (2019) showed that vertical eddy diffusion coefficients due to small-scale turbulences should be less than 0.02 m² s⁻¹ to reproduce realistic super-rotation. This value is significantly smaller than the eddy diffusion derived using the above formulation. This could suggest that global-scale circulation effectively represents the eddy diffusion transport in a 1D column. Specifically, the dynamical timescales obtained in the GCMs are consistent with the timescale required to reproduce observed cloud structure and condensational gas species. Above 85 km, the existence of subsolar-to-antisolar (SSAS) circulation is supported by observations of O₂ nightglow at the antisolar point (Soret et al. 2012) and the Venus Thermospheric GCM (e.g., Brecht et al. 2012; Hoshino et al. 2013; Navarro et al. 2021). Thus, the large eddy diffusion above 85 km may be attributed not only to Hadley-type circulation but also to the influence of the SSAS circulation.

The case studies showed that the H₂O VMR profile is highly sensitive to the eddy diffusion coefficient between 60 and 70 km. The eddy diffusion coefficient used in this study is broadly consistent with previous observations and numerical models. Therefore, a thorough examination of the eddy diffusion coefficient in this altitude range is crucial for understanding H₂O transport through the cloud layer. Numerical investigations of the eddy diffusion coefficient can be performed using a GCM, as shown by Zhang & Showman (2018). Recently, Fujisawa et al. (2022) achieved a realistic representation of wind and thermal structures using a data assimilation technique that combines AFES-Venus and AKATSUKI observations. Such data assimilation products could serve as valuable resources for deriving an eddy diffusion coefficient that accounts for both global and local transport processes.

In terms of the mesospheric temperature, the SOIR profile is more suitable to reproduce the observed H₂O VMR increase above 85 km (see Figure 13). Therefore, the present study suggests that Case 5 provides the best match with observed atmospheric parameters. However, observed temperature profiles appear to differ depending on the instrument used for observation. In general, radio occultation measurements yield lower temperature profiles compared to those obtained through infrared observations (Limaye et al. 2018), as depicted in Figure 3(a). In addition, Krasnopolsky (2011) argued that high temperatures between 95 and 100 km disagree with the rotational temperature of the O₂ nightglow at 1.27 μm.

This discrepancy could potentially be resolved by measuring the liquid H₂SO₄ mass fraction in upper haze particles. As shown in Figure 9(a), the liquid H₂SO₄ mass fraction varies significantly depending on the temperature—a variation that should be detectable. Consequently, we can use the H₂SO₄ mass fraction as a proxy to better understand the upper-mesospheric temperature profiles and the resulting equilibrium conditions. However, as Mahieux et al. (2023b) pointed out, the temperature difference could be caused by the local time dependence. In that case, the temperature observations can be used as a proxy for the H₂SO₄ mass fraction of aerosols. Regardless of the scenario, acquiring more detailed observations that focus on the local time dependence of mesospheric temperatures is crucial for enhancing our understanding of the microphysical properties of aerosols.

The H₂O abundance at the homopause regulates the hydrogen (H) escape rate through diffusion-limited flux. Assuming a H₂O vapor abundance of ∼1 ppm and HCl of ∼0.5 ppm at the homopause, the diffusion-limited flux is calculated to be ∼3 × 10¹¹ m⁻² s⁻¹ (Catling & Kasting 2017). This value is consistent with an estimate of the nonthermal escape flux of ∼3 × 10¹¹ m⁻² s⁻¹ (Lammer et al. 2006), indicating that the diffusion-limited flux serves as a good approximation for estimating the present-day H escape flux. Therefore, the H₂O vapor abundance at the homopause is crucial for understanding the current atmospheric evolution of Venus.

The long-term stability of the Venusian clouds is still controversial. Bullock & Grinspoon (2001) showed that the clouds would not be stable if the atmospheric SO₂ amount is buffered by carbonate. Conversely, Hashimoto & Abe (2005) argued that the clouds would be stable if the atmospheric SO₂ amount is buffered by pyrite. These contrasting scenarios could lead to significant differences in H₂O loss rates, given that the cloud layer affects the diffusion-limited escape rate. Therefore, the presence of clouds and their impact on mesospheric H₂O abundance should be carefully considered in models of Venusian atmospheric evolution. The present study shows that cloud processes are essential to evaluate the H₂O VMR vertical profile, which is deeply connected to H₂O vapor abundance at the homopause. Consequently, cloud processes cannot be disregarded when studying the mesospheric H₂O abundance and the diffusion-limited escape rate. These cloud effects should be considered not only for the present-day atmosphere but also for the past atmospheric conditions. Furthermore, observations of trace gas species and clouds, possible with the future EnVision mission, will further elucidate the effects of clouds on atmospheric evolution (e.g., Widemann et al. 2023). The VenSpec-H instrument on board the Envision spacecraft will observe H₂O abundance below and above the cloud layer with high accuracy, and the Radio Science experiment will obtain liquid H₂SO₄ abundance in the cloud layer (Akins et al. 2023). This information will allow the community to investigate the role of clouds in regulating the mesospheric H₂O abundance and its impact on the hydrogen escape.

The simulated size distribution in the upper haze layer is unimodal, diverging from the bimodal distribution suggested by VEx observations (Wilquet et al. 2009). Previous microphysics studies (Gao et al. 2014; Parkinson et al. 2015) have also reported unimodal distributions in the upper haze layer at steady state. However, Gao et al. (2014) argued that transient vertical winds could give rise to the observed bimodal distribution. The inclusion of dynamical effects, which is not considered in the present study, may be necessary to accurately reproduce the observed particle size distribution and should be explored in future research.

Homogeneous nucleation processes could play a significant role in the generation of the Mode 1 particles. Dai et al. (2022) obtained a supersaturated H₂SO₄ VMR profile in the upper cloud and haze layer, which aligns with the findings of our study; their results were obtained using a cloud condensation model. They argued that under such a supersaturated condition homogeneous nucleation is likely to occur. The enhanced CN population could correspond to Mode 1 particles, as observed in the lower and middle cloud layers of the recent microphysical model (McGouldrick & Barth 2023). This idea is supported by the VIRA temperature cases (Nominal Case and Cases 1–3), where the H₂SO₄ vapor is highly supersaturated in the upper haze layer, as shown in Figure 7(a). The potential for homogeneous nucleation should be further explored using a cloud microphysics model that includes CN size distributions in future studies.