The Impact of Turbulent Vertical Mixing in the Venus Clouds on Chemical Tracers

Venus clouds host a convective layer between roughly 50 and 60 km that mixes heat, momentum, and chemical species. Observations and numerical modelling have helped to understand the complexity of this region. However, the impact on chemistry is still not known. Here, we use for the first time a three-dimensional convection-resolving model with passive tracers to mimic SO$_2$ and H$_2$O for two latitudinal cases. The tracers are relaxed towards a vertical profile in agreement with measured values, with a timescale varying over several orders of magnitude. The vertical mixing is quantified, it is strong for a relaxation timescale high in front of the convective timescale, around 4 hours. The spatial and temporal variability of the tracer due to the convective activity is estimated, with horizontal structures of several kilometres. At the Equator, the model is resolving a convective layer at the cloud top (70 km) suggested by some observations, the impact of such turbulent activity on chemical species is accounted for the first time. From the resolved convective plumes, a vertical eddy diffusion is estimated, consistent with past estimations from in-situ measurements, but several orders of magnitude higher than values used in 1D chemistry modelling. The results are compared to observations, with some spatial and temporal variability correlation, suggesting an impact of the convective layer on the chemical species.


Introduction
The strong dynamical activity inside the Venusian cloud layer has been assessed since the beginning of the Venus spacecraft exploration. Cloud images by the Mariner 10 mission (Belton et al., 1976) and the Pioneer Venus spacecraft (Rossow et al., 1980) near the subsolar point showed cellular features suggesting convective cells with diameters between 200 and 1000 km. The convective activity was first measured by the Pioneer Venus radio occultation experiment (Seiff et al., 1980) from 50 to 55 km of altitude, and was then confirmed by the radio occultation experiment onboard the Magellan probe (Hinson and Jenkins, 1995). The VeGa balloons flew in that altitude range close to the Equator and measured vertical winds between -4 and 2 m s´1 Lorenz et al., 2018) and convective cell diameter from several hundred meters to tens of kilometers  around 54 km of altitude. The VeRa radio occultation device on board of Venus Express studied in detail this convective layer and measured a strong latitudinal variability of the depth of the layer (Tellmann et al., 2009), reaching 10 km close to 80˝of latitude, almost twice the value of the equatorial regions. The radio occultation experiment in the Akatsuki spacecraft measured variability of the convection depth with local time , this layer being thicker a night.
In addition to the convection layer in the deep cloud layer, the Venus Monitoring Camera (VMC) observed the cellular features at the top of the cloud, about 70 km of altitude, at the subsolar point suggesting convective activity (Markiewicz et al., 2007;Titov et al., 2012). A convective layer at this altitude is the main hypothesis for these observed structures, with measured convective cells from 20 to a few hundred of kilometres. However, the different radio occultation on board of Venus Express and Akatsuki radio occultation did not measure any clear neutral-stability layers at the subsolar point (Ando et al., , 2020. Gravity waves emitted from the convective layer have been observed by different space probes, Pioneer Venus radio science observed evidenced small-scale waves with vertical wavelengths of about 7 km above and below the cloud layer (Seiff et al., 1980;Counselman et al., 1980), the Venus Express instruments measured the wavelengths of the waves emitted above the cloud layer, ranging between about 2 and 3.5 km vertically (Tellmann et al., 2012) and from 2 km to hundreds of kilometres horizontally (Peralta et al., 2008;Piccialli et al., 2014). The gravity waves in this region have also been studied with Akatsuki Mori et al., 2021).
Decades of spacecraft and ground-based observations of sulphur dioxide and water show highly variable abundance in the upper cloud deck, with timescales from hours to decades Encrenaz et al., 2016;Vandaele et al., 2017a,b). Convection is one of the hypotheses for the short, from hours to days, term variability (Marcq et al., 2013;Vandaele et al., 2017a). HST Imaging Spectrograph was used to observe Venus (Jessup et al., 2015 at cloud-top altitudes, albedo darkening was measured and explained by a possible increase of the convective vertical mixing and the injection of the unknown absorbing species. Using the SOIR/Venus Express CO 2 and CO profiles above the clouds and 1D photochemical model, Mahieux et al. (2021) estimated the vertical mixing from 80 to 140 km.
1D models have been developed to study the chemistry of the Venusian atmosphere. Krasnopolsky (2007) and Krasnopolsky (2013) focused on the lower atmosphere and Krasnopolsky (2012), Zhang et al. (2012) and Parkinson et al. (2015) and Shao et al. (2020) on the middle atmosphere. Yung et al. (2009) and Bierson and Zhang (2020) and Rimmer et al. (2021) modelled the atmosphere from the surface to 110 km. These models cannot resolve the turbulent activity inside the cloud and therefore use the eddy diffusivity coefficient formalism to represent the different turbulent processes in the atmosphere like convection. There is a large uncertainty over the value of this coefficient in the Venusian atmosphere. These models showed that the mesospheric abundance of several species (especially SO 2 ) was very sensitive to the eddy diffusivity values and vertical gradient.
Due to the lack of understanding of the turbulence inside the clouds of Venus, the effect of the cloud convective layer and gravity waves on the chemistry and microphysics has not been studied in detail. Only McGouldrick and Toon (2008) gave an insight into the change of optical depth due to the convection and gravity waves, using an idealized 2D (zonal/vertical) representation of the Venus cloud convective layer. Morellina and Bellan (2022) studied the vertical mixing due to the species stratification in the Venus lower atmosphere and clouds using Direct Numerical Simulation. High density-gradient magnitude regions are formed with increasing stratification and low stratification conditions produce a more uniform spatial distribution of the density.
To understand the turbulence activity inside the Venus cloud layer, the limited-area Venus mesoscale model (VMM) adapted from a terrestrial hydrodynamical solver (Skamarock and Klemp, 2008) was developed at Laboratoire de Météorologie Dynamique by Lefèvre et al. (2017) and later coupled to the full set of physical packages for Venus developed at Institut Pierre Simon Laplace (IPSL) (Lebonnois et al., 2010(Lebonnois et al., , 2016Garate-Lopez and Lebonnois, 2018) to simulate only a specified region of the planet with a fine resolution.
The 3D Large-Eddy Simulation (LES) mode was used to study the convection and smallscale gravity waves at different latitudes and local times (Lefèvre et al., 2017(Lefèvre et al., , 2018. The resolved convection depth is consistent with observations with updrafts cells diameter of 20 km. The gravity waves emitted by the convective region are also consistent in amplitude and wavelength with observations. Cloud-top convective activity is also present around the subsolar point.
In this study we propose to use this convection-resolving to model at noon for two latitudinal cases, the Equator and 75˝, with idealized passive tracers representing SO 2 and H 2 O to quantify the impact of the resolved convective layer vertical mixing on the Venusian cloud layer chemistry, and to compare the results with 1D chemical models and observations. The resolved convective plumes give an unprecedented insight into the cloud convective advection and chemistry dynamics.
Our paper is organized as follows. The model is described in Section 2. In Section 3, the impacts of convective motions and gravity waves are presented. The vertical eddy diffusivity is estimated in Section 4. Results are discussed in Section 5. Our conclusions are summarized in Section 6.

Mesoscale modelling for Venus
Our mesoscale model for Venus is based on the dynamical core of the Advanced Research Weather-Weather Research and Forecast (hereinafter referred to as WRF) terrestrial model (Skamarock and Klemp, 2008). The WRF dynamical core integrates the fully compressible non-hydrostatic Navier-Stokes equations on a defined area of the planet. The conservation of the mass, momentum, and entropy are ensured by an explicitly conservative flux-form formulation of the fundamental equations (Skamarock and Klemp, 2008), based on mass-coupled meteorological variables (winds and potential temperature). The LES mode is used in this study, and the parametrization of the unresolved small-scale eddies is performed using a subgrid-scale prognostic Turbulent Kinetic Energy closure by Deardorff (1972). This method has been used for Earth convection study (Moeng et al., 2007), for the Martian atmosphere (Spiga et al., 2010), and for terrestrial exoplanets .
Following the work of Lefèvre et al. (2018) for Venus, the WRF core is coupled to the radiative transfer code of the IPSL Venus General Circulation Model (GCM) (Lebonnois et al., 2015) based on Eymet et al. (2009) net-exchange rate (NER) formalism computing the energy between the layers prior to the dynamical simulations, separating temperature-independent coefficients from the temperature-dependent Planck functions of the different layers. The cloud model is based on Haus et al. (2014) and Haus et al. (2015), with a latitudinal variation of the cloud by setting 5 distinct latitude intervals: 0˝to 50˝, 50˝to 60˝, 60˝to 70˝, 70˝to 80a nd 80˝to 90˝. The cloud model is set prior to the simulations and does not interact with the dynamical features resolved by the model. The solar heating rates are computed using look-up tables (Haus et al., 2015) of vertical profiles of the solar heating rate as a function of the solar zenith angle.

A Passive tracer approach
Tracers have been included in the model and are set to represent SO 2 and H 2 O . The tracers are not radiatively active. The chemistry, photodissociation, and condensation sources and sinks are modelled by a linear relaxation of the tracer abundance q toward a prescribed vertical profile q 0 pzq with a relaxation timescale τ computed as : The two prescribed relaxation profiles are displayed in Fig 1. These two profiles are constructed as follows: assuming a constant abundance value below the clouds from spectroscopic measurements between 30 and 40 km (Bézard and de Bergh, 2007), 150 ppm for SO 2 and 30 ppm for H 2 O . This value of SO 2 is consistent with in-situ measurements at the bottom of the cloud deck by VeGa 1 and VeGa 2 entry probes (Bertaux et al., 1996). From 48 km upwards an exponential decay is assumed going through 0.5 ppm at 65 km, long-term ground observations , for SO 2 and 2 ppm at 70 km for H 2 O , VIRTIS and SPICAV measurements (Cottini et al., 2012;Fedorova et al., 2016). These profiles have been constructed to represent values in the cloud layer, and were constructed for simplicity as constant below the cloud deck and decreasing above. With such simple profiles, there are discrepancies with observations. The value of SO 2 in the lower cloud is slightly underestimated (Oschlisniok et al., 2021). Above the clouds, both SO 2 and H 2 O are decreasing despite the inversion observed in SPICAV/SOIR  and SPICAV/UV (Evdokimova et al., 2021) for SO 2 , and the very low vertical gradient between 80 and 120 km for H 2 O observed by SOIR/VEx (Chamberlain et al., 2020). The profile of both SO 2 and H 2 O outside the cloud layers does not impact the results obtained within the clouds. No source or sinks at the surface or the top boundary are considered. The lateral boundary conditions are doubly periodic. Two latitudinal cases are considered in this study, the Equator and 75˝, to explore the latitudinal variability in terms of convective activity. The relaxation profiles are the same for the two latitudinal cases. Only one local time is considered in this study because of the small diurnal cycle of the convective layer in the model. At night, there is no cloud-top convective layer and this region is similar in terms of dynamics to the high latitude at noon. The chemical timescale of SO 2 and H 2 O is not well-constrain over the column, latitudes and local time. Therefore, the relaxation timescale τ is set to a constant value over the column, with sensitivity tests ranging over values of 10 2 , 10 3 , 10 4 , 10 5 and 10 6 s. This broad range of relaxation timescale was chosen to be compared to the dynamical timescale of the convective layer around 10 4 s (see Section 2.3) with two orders of magnitude above and below this value. Shao et al. (2020) reports chemical timescale below 10 4 s in the upper cloud and as high as 10 8 s. Zhang et al. (2012) reports that the nucleation timescale in the clouds can be below 10 2 s. There are no bottom or top boundary conditions for the tracers, the initial profiles for SO 2 and H 2 O are considered at equilibrium in regard to chemistry, photochemistry and condensation/vaporization.

Simulation settings
The simulation settings in terms of resolution and time-step are identical to Lefèvre et al.  Fedorova et al. (2016), the star represents the SO 2 value from Encrenaz et al. ( , 2015 and the horizontal line represents tropospheric abundance measurements from Bézard and de Bergh (2007).
with in-situ measurements of the VeGa balloons Sagdeev et al., 1986). Regarding the gravity waves, with or without the presence of a cloud-top activity, the amplitude of the gravity wave in temperature perturbations and vertical wavelengths are consistent with the radio-occultation measurements (Tellmann et al., 2012), and the horizontal wavelength is consistent with cloud-top UV observations (Piccialli et al., 2014). With a realistic radiative transfer and incoming solar heating, characteristic of an average UV absorber abundance , the model exhibits a cloud-top convective layer with convective cells diameters consistent with the VMC observations (Markiewicz et al., 2007;Titov et al., 2012) that is still speculative but the main hypothesis due to the absence of measurements of corresponding static stability vertical distribution. This cloud-top convective layer has only a limited impact on the deep cloud convective layer because the source, IR heating at cloud base from the troposphere, is unaffected. The thin neutral at 75°of latitude at 7 km is an artefact from the radiative transfer and large-scale heating extracted from the LMD Venus GCM with a much more coarse resolution than the convection-resolving model. To avoid spurious reflection of upward propagating gravity waves on the top boundary of the model, a Rayleigh damping layer is applied over the 8 top kilometres with a damping coefficient of 0.08 s´1. The tracers are initialized with relaxation profiles (Fig 1) and then advected by the dynamics: the convection, wind shear, and gravity waves. The outputs of the simulations are shown after about two Earth days of simulation.

Mixing
Fig 3 shows the domain averaged vertical profiles of SO 2 at the Equator and 75˝. The presence of the convective layer, well-mixed, in the deep cloud region is noticeable with a constant tracer abundance value for larger relaxation timescales of 10 5 and 10 6 s (Fig 3-b). To illustrate the vertical mixing in the cloud layer, a dynamical timescale is defined as follows τ dyn = H/σ w , with H being the scale height and σ w the spatial standard deviation of the vertical wind at a given altitude, and τ dyn represents the mixing timescale. At the Equator, τ dyn is equal to about 17000 s at 50 km and equal 6000 s at 70 km. At high latitude τ dyn is equal to 8000 s at 50 km.
For a relaxation timescale of 10 4 s, the value of the tracer is not constant in the convective layer but has a significantly steeper vertical gradient compared to the initial state. There is hardly any discernible difference between the initial profile and the tracer abundance profiles for relaxation timescale of 10 2 and 10 3 s. When the relaxation timescale is significantly larger than the dynamical timescale τ dyn , like for the 10 5 and 10 6 s cases at the Equator, the changes in abundance due to the convection acting faster than the chemical relaxation, and therefore the tracer will be well-mixed. The term chemical relaxation encompasses actual chemical reactions, but also condensation and sublimation for H 2 O .
On the contrary, when the relaxation timescale is smaller by several orders of magnitude than the dynamical timescale τ dyn , in the 10 2 and 10 3 s cases for example, the chemistry will act faster than the perturbations in abundance due to the convection and the abundance will tend to the prescribed value. When the dynamical timescale τ dyn and the relaxation timescale are close, like in the 10 4 s case, the tracer abundance will be in-between the two extreme cases. The same competition between the chemical/physical equilibrium and the well-mixed case is also visible at the cloud top. At the Equator, the cases with a relaxation timescale greater than the 6000 s dynamical timescale τ dyn are well-mixed.
At 75˝of latitude the deep convective layer is slightly thicker with stronger vertical wind, therefore the dynamical timescale τ dyn is smaller than at the Equator, thus there is a slight difference for the in-between 10 4 s case, the tracer being closer to the well-mixed case at high latitude. At the cloud top, there is no convective activity, the dynamical timescale becomes very large, and the tracer is always close to the prescribed profile.
Other noticeable features are the abundance drop at the top of the two convective layers and the smaller increase at the bottom of the deep convective layer. These two features are due to the convective entrainment layers, dominated by the updrafts above the convective layers and by the downdrafts below (See Fig 6 of Lefèvre et al. (2017)). The amplitude of these features increases with increasing relaxation timescale; when the convection dominates the vertical mixing, the convective overshoots will have a greater impact.

Spatial variability
Fig 4 displays snapshots of SO 2 abundance at the same three altitudes in the cloud and for four cases of relaxation timescale at the Equator. The 10 6 s case is not shown here for clarity, since it yields very similar results to the 10 5 s case. The spatial variability is smaller for extreme relaxation timescale, 10 2 and 10 5 , resulting in smaller intervals ranging less than 10% in relative value. At 52 km (left column), the effect of the upward convective plume is visible by the high abundance values forming structures up to 10 km in diameter and narrow in-between. At 62 km (middle column), the differences between the different relaxation timescales are hardly noticeable. However, the gravity wave horizontal structure is well-perceptible. As for the cloudtop convective layer (right column), the relaxation timescale will impact in the same way the interval values and the horizontal structure. To precisely measure the horizontal structure of chemical species in the clouds, a resolution of the order of the kilometre is necessary. These trends are similar for H 2 O (see Fig 9 in Appendix B), and at 75˝as well, with the presence of gravity waves instead of a convective layer at 70 km for SO 2 (Fig 8 in Appendix A).   at three altitudes, 52, 62 and 70 km representing respectively the deep convective layer, the gravity waves regions and the cloud-top altitudes, and for both species SO 2 (left) and H 2 O (right). The competition between the chemical/physical equilibrium and the convective mixing is readily visible in this figure. At 52 km, As shown in Fig 4, for the extreme relaxation timescale cases, 10 2 , 10 5 and 10 6 s, the predominance of either the chemistry or the convective mixing over the abundance leads to a spatial relative standard deviation below 0.1, lower by a factor 10 than the in-between cases. Any spatial perturbation would be either mixed by the convection or converted by the chemistry, through a sink or a source. At 62 km, the vertical mixing ensured by the gravity waves is comparatively small (see Fig 2), and therefore the spatial abundance disparities are smaller. At 70 km, when there are gravity waves the spatial relative standard deviation is small and when there is the cloud-top convective layer, the spatial relative standard deviation will be maximum for the in-between cases like at 52 km. The relaxation timescale also has an impact on the vertical gradient of the tracers inside the convective layer. For the extreme relaxation timescale cases, 10 2 , 10 5 and 10 6 s, the vertical gradient will either be strong when the chemistry dominates, or be almost zero when the convective mixing dominates. For intermediate relaxation timescale, there can be the presence of non-monotinic vertical gradient inside the convective layer.

Temporal variability
Fig 6 shows the temporal relative standard deviation for the latitudinal cases for the same three altitudes, 52, 62 and 70 km for SO 2 (left) and H 2 O (right). The time step is determined as minpτ dyn , τ q{10, and the standard deviation is computed over a hundred time steps. The competition between the chemical/physical equilibrium and the convective mixing is visible and is of the same order as that of the spatial variability (Fig 5). The temporal variability is higher, superior to 0.1, for the 10 3 and 10 4 s cases in the convective layers, at least one order of magnitude over the other location. The physical interpretation is the same as for the spatial variability, related to the competition between τ dyn and τ .

Vertical eddy diffusion
The vertical eddy diffusivity, or K zz , is representing the vertical mixing by turbulent processes such as convection or gravity waves, and is used in atmospheric modelling as a parametrization when the turbulence is not resolved, like in 3D GCMs or 1D modelling. For LES models, like the one presented here, with a small resolution, the larger eddies are resolved and no vertical eddy diffusivity is therefore needed.
K zz has been estimated in the Venus atmosphere from observations and modelling, but is still not well-constrained, ranging over several orders of magnitude. The Pioneer Venus radio scintillation measurements estimated K zz between 2¨10´1 m 2 s´1 at 45 km (Woo et al., 1982) and 4 m 2 s´1 at 60 km (Woo and Ishimaru, 1981). Using the values of the vertical wind measured in the deep cloud layer by the VeGa balloons, of the order of one meter per second , K zz was evaluated at 10 3 m 2 s´1 at 54 km  and between 1 and 10 2 m 2 s´1 in the 50-56 km range (Imamura and Hashimoto, 2001;Gao et al., 2014). McGouldrick and Toon (2007) estimated K zz between 1 and 10 3 m 2 s´1 in the 50-57 km range. Above the clouds, K zz was estimated using the SOIR/Venus Express CO 2 and CO profiles at around 360 m 2 s´1 from 80 km (Mahieux et al., 2021).
As the vertical plumes are resolved in LES simulations such as ours, K zz can be estimated from tracers' variations as follows: with q a tracer mixing ratio, primed quantities representing perturbations relative to the domain averaged values, and bracket quantities representing domain averaged values. Fig 7 displays the vertical eddy diffusivity profiles for the two latitudinal cases and the four relaxation timescale considered with SO 2 . Below the clouds, the value of K zz is in the order of 10´1 m 2 s´1 consistent with (Woo et al., 1982). In the deep convective layer, K zz increases and reaches values up to 10 3 m 2 s´1 at the Equator and 10 4 m 2 s´1 at high latitudes. The order of magnitude is consistent with estimation using Prandtl mixing-length theory (Lindzen, 1971). The values for relaxation timescale equal to 10 4 , 10 5 and 10 6 s at the Equator are consistent with in-situ measurements from VeGa balloons at the same latitudes . Most of the values between 50 and 57 km are consistent with the estimation of Imamura and Hashimoto (2001); Gao et al. (2014) and McGouldrick and Toon (2007). However, at high latitudes the vertical eddy diffusivity values exceed these estimations, reaching 10 4 m 2 s´1, due to a thicker convective layer, consistent with VeRa measurements (Tellmann et al., 2009). The relaxation has an important impact on the vertical mixing, and changes the value of the vertical eddy diffusivity by several orders of magnitude, for high τ values the tracer is wellmixed with a small vertical gradient in the convective layer, resulting in a high vertical eddy K zz (m 2 s´1)

Previous estimations
This study at 45 km 2¨10´1 (Woo et al., 1982) 10´1-1 at 54 km 10 3  50-4000 50-56 km 1-10 2 (Imamura and Hashimoto, 2001) 20-10 000 50-57 km 1-10 3 (McGouldrick and Toon, 2007) 10´1-10 000 at 60 km 4 (Woo and Ishimaru, 1981) 10´1-20 diffusivity. Above, K zz decreases to 1 m 2 s´1, consistent with in-situ measurements (Woo and Ishimaru, 1981) at 60 km. At noon near the Equator, the gravity waves resolved by the model are trapped between the two convective layers and have a stronger amplitude. The vertical eddy diffusivity is therefore stronger than at higher latitudes, with an opposite effect upon the relaxation timescale, where K zz is stronger for a lower value of τ . At the cloud top at noon, the vertical eddy diffusivity strongly increases to values up to 10 4 m 2 s´1 due to the cloud-top convective layer generated by solar heating absorption from the unknown UV absorber (Lefèvre et al., 2018), where the relaxation has an important effect similarly as the deep cloud convective layer. The presence of a cloud-top convective layer is still speculative, however the turbulence resolved by the model provides an insight of the impact of such convective activity and trapped gravity waves. Without the presence of this cloud-top convective layer, the mixing of the gravity waves at the Equator above the deep convection would be equivalent to the mixing at 75˝of latitude in the gravity waves region. Regarding the local time variability, the deep convective layer is deeper at night Imamura et al. (2017); Lefèvre et al. (2018). Therefore, the mixing would be stronger. The small increase for 75˝of latitude at 72 km is an artefact mentioned in Section 2.3. The vertical eddy diffusivity profiles were also estimated using H 2 O , and the orders of magnitude are consistent between the two species. The estimation of the vertical eddy diffusion for heat, using potential temperature instead of tracer abundance in Equation 2, gives a value between 10 3 and 10 4 m 2 s´1 for the Equator in the deep cloud convective region and around 10 4 m 2 s´1 for the cloud-top convective layer, and between 10 4 and 10 5 m 2 s´1 for high latitudes. Above the clouds and below the Rayleigh damping layer, K zz is estimated from 1 to 10 m 2 s´1, at least one order of magnitude lower than the values estimated from 80 km (Mahieux et al., 2021). The small-scale turbulence only plays a minor role in the mixing in this region. The present estimation of the vertical eddy diffusivity coefficient and comparison with previous estimations are summarized in Table 1. Figure 7: Vertical profiles of the vertical eddy diffusivity (m 2 s´1) in the Venus cloud region calculated with the SO 2 tracer for relaxation timescale of 10 2 s (green), 10 3 s (blue), 10 4 s (yellow), 10 5 s (cyan) and 10 6 s (red) for the Equator (solid line) and 75˝(dashed line). The black circles represent the previous estimations (see Table 1) and the black line represents the range of values used for the cloud convective layer in recent chemical models Krasnopolsky (2012); Bierson and Zhang (2020); Rimmer et al. (2021).

Discussion
The Krasnopolsky (2012) 1D model uses an eddy diffusion of 1 m 2 s´1 in the convective layer, at least 10 times lower than the present study, constant up to a specific altitude set to 55, 60 or 65 km, and then increases linearly. Both modelled SO 2 and H 2 O abundance are sensitive to small variations of eddy diffusion, not constant in the convective layer, suggesting a short chemical timescale, inferior to 10 4 s. The 55 and 60 km top altitude of the deep eddy diffusion constant value results in a better agreement with mesospheric measurement values. This eddies diffusion constant value depth is consistent with the convective layer depth in this study. The 1D model of Bierson and Zhang (2020) defined several eddy diffusion scenarios for the convective layer mixing parametrization, from 0.1 to 2 m 2 s´1 nominal eddy diffusion value, at least one order of magnitude lower than the present study. With the nominal eddy diffusion set-up, constant from the surface to 60 km, the transport timescale ranges from 10 years at the surface to a few months at 90 km, orders of magnitude larger in the convective layer than the convective dynamical timescale τ dyn . However, SO 2 exhibits in their model a very well-mixed profile below 60 km, a step decrease and then another region with a very low vertical gradient up to 110 km, inconsistent with the observations from Venus Express Vandaele et al., 2017a;Evdokimova et al., 2021). The exploration of the eddy diffusion values and profile in the convective layer leads to a better agreement with the observations of SO 2 abundance in the mesosphere with a lower value. Rimmer et al. (2021) also tested several vertical eddy diffusion values in the convective layer, between, used also 0.01 to 1 m 2 s´1. However, these tested eddy diffusion values are inconsistent with estimations from in-situ wind and cloud microphysics, and the values presented in Section 4.
The chemical models of Yung et al. (2009) andZhang et al. (2012) and Shao et al. (2020) simulate the atmosphere from 58 km, therefore above the convective layer. Shao et al. (2020) showed that the temporal variations of SO 2 and H 2 O are linked to the lower atmospheric processes, such as convection. No chemistry model considers the cloud-top convection hypothesis through the inclusion of variable K zz profiles. In the IPSL Venus GCM, K zz reaches values above 10 m 2 s´1 between 48 and 57 km (Lebonnois and Schubert, 2017), the lowest value estimated with the present convection resolving model.
The minor gaseous species have been measured mainly in the mesosphere by Venus Express above 60 km, e.g. Marcq et al. (2020) for SO 2 and Cottini et al. (2012);Fedorova et al. (2016); Vandaele et al. (2017a); Chamberlain et al. (2020) for H 2 O , and around 60 km by ground-based telescopes (Encrenaz et al., 2019 and Hubble Space Telescope (Jessup et al., 2015. The spatial resolution of Venus Express VIRTIS and SPICAV instruments can be around a few tens of kilometre at best (Drossart et al., 2007;Bertaux et al., 2007), while the spatial resolution is about 100 km for the ground-based observations and between 20 and 60 km for the Hubble measurements. The uncertainty about the height of the probed regions for both IR and UV observations makes the comparison between the present model and observations delicate.
A short-time variability of two hours for SO 2 was measured from the ground (Encrenaz et al., 2016) at around 60 km, consistent with the dynamical timescale of the convective layer presented here. The probability of this variability is stronger at night (Encrenaz et al., 2019), consistent with the local variability of the convective layer in the model (Lefèvre et al., 2018) and observations . VeGa 1 and VeGa 2 entry probes measured SO 2 between 60 km and the ground (Bertaux et al., 1996). Inside the convective layer, a strong nonmonotonic vertical gradient was measured, consistent with the vertical gradient for intermediate relaxation timescale of SO 2 in the same region but with a higher amplitude.
At cloud-top altitudes, observed SO 2 does not show clear signs of a vertically well-mixed layer with a large set of values of the vertical gradient (Vandaele et al., 2017a), suggesting a short chemical timescale (inferior to 10 3 s) or no cloud-top convection layer. However, SO 2 was measured at the top of the cloud layer by SPICAV-UV, a higher variability is observed below 30˝of latitude with an order of magnitude consistent with the presence of a cloud-top convective layer , in the same latitude range as the cellular features observed by VMC (Titov et al., 2012), consistent with a relaxation timescale between 10 3 and 10 4 s. At these latitudes, the abundance of observed SO 2 is smaller around noon correlated with an increase of the unknown UV absorber, whose solar heating is needed to generate the cloud-top convective layer. This minimum of SO 2 at noon was also observed from the ground (Encrenaz et al., 2019. In the same altitude range, H 2 O was measured by VEX/VIRTIS in the turbulent-like cloud features (Cottini et al., 2012). The horizontal abundance variation in this area is about 1 ppm, consistent with our results (Fig 5 and Fig 9) with a relaxation timescale of 10 3 and 10 4 s.
SPICAV measured H 2 O in the upper cloud (Fedorova et al., 2016), between 61 and 59 km (the lowest altitude probed) the vertical gradient is very weak, consistent with a well-mixed layer. A small vertical gradient of H 2 O is needed for better retrievals of HDO ground observations (Encrenaz et al., 2013), whereas a significant vertical gradient is used for SO 2 . (Encrenaz et al., 2016) measured a higher variability for SO 2 than for H 2 O , suggesting again a well-mixed layer for water. The variations observed at the cloud top show no signs of a diurnal cycle at low latitudes, suggesting also a short chemical timescale (inferior to 10 3 s) or no convection layer.

Conclusion
This study is the first one to investigate the vertical mixing in the Venusian cloud layer on the chemistry using 3D resolved convective plumes. Using an idealized passive tracer to represent SO 2 and H 2 O , vertical mixing of those species by the convection was estimated depending on the chemical timescale. A small chemical timescale compared to the convective timescale, about 10 3 s or lower, will limit the vertical turbulent mixing. On the contrary, a chemical timescale value superior to the convective timescale, about 10 5 s or higher, will lead to a wellmixed layer, i.e. a small vertical gradient in the convective region. The convective region is organized as updraft polygonal cells, and the spatial and temporal variability through the cloud is estimated. These variabilities are maximal in the convective regions for a convective chemical timescale equivalent to the convective timescale. The horizontal structure of the SO 2 and H 2 O is composed of small-scale structures due to the convective cells or gravity waves as small as a few kilometres in the convective regions. In the gravity wave area, the horizontal structure of the tracer is governed by the gravity wave horizontal wavelength. The mixing on tracers of the hypothetical cloud-top convective layer has been estimated, with a smaller convective timescale, the convective plumes can engender small-scale structures. With the presence of such a cloudtop convective layer, the gravity waves have larger amplitudes and exhibit a stronger mixing.
The higher values of the mixing in the upper cloud layer could help to infer the presence of a cloud-top convective layer.
From the resolving plumes, the vertical eddy diffusion has been estimated in the cloud region. For the deep cloud convection, K zz ranges from 10 1 to 10 4 m 2 s´1, consistent with in-situ and modelling estimation. There is an impact of the relaxation timescale, the stronger it is, the higher K zz is, and an impact of the latitudinal convection depth variability. In the gravity waves region, the vertical eddy diffusion is lower, and increases again at the cloud top at the subsolar point due to the cloud-top convective activity.
Previous 1D chemical models Yung et al. (2009); Bierson and Zhang (2020); Rimmer et al. (2021) found no mixed layers for SO 2 and H 2 O . However, the various vertical eddy diffusion profiles and values tested are not consistent with the values of the present study and convective mixing is underestimated. These models show that vertical eddy diffusion is important for the mesospheric abundance of the chemical species, there is a need for new 1D chemical simulations with a more realistic convective mixing.
The limited size of the domain and the simplicity of the representation of the chemistry makes comparison with the observations difficult. Nevertheless, there are interesting points to comment on. A two-hour variability of SO 2 around 60 km measured from the ground is compatible with the dynamical timescale of the convective layer. The horizontal abundance variation of H 2 O in the turbulent-looking cloud top is consistent with the spatial variability of the subsolar cloud-top convection. The weak vertical gradient of H 2 O between 61 and 58 km is consistent with convective layer vertical mixing. The vertical profiles of SO 2 and H 2 O do not indicate clear signs of cloud-top convective mixing. However, the cloud-top variability of SO 2 around the subsolar point is consistent with the presence of a cloud-top convective layer.

Future improvements
A wider horizontal domain is needed to be able to compare observations with more accuracy. The passive tracer method used in this study is idealized, a proper chemistry and microphysics scheme are needed for a more realistic model, and to be able to quantify the impact of the convective activity on cloud formation and cloud opacity or to study the impact of turbulence on coupled particles trend, like the anti-correlated variations of H 2 O and SO 2 predicted by 1D chemical model (Parkinson et al., 2015;Shao et al., 2020) and observed from the ground .
The cloud-top convective activity resolved in the model is generated by the unknown UV absorber. However, its abundance distribution is not well-constrained and this uncertainty could have a strong effect on the cloud-top convection depth and vertical mixing. The observed variations of this unknown absorber  need to be taken into account for a better understanding of this cloud-top convective layer.
This study focuses on small-scale turbulence, a similar study needs to be conducted for the mesoscale features observed in the clouds with Akatsuki, like the vertical propagation mountain (Kouyama et al., 2017;Lefèvre et al., 2020) and the other new cloud morphologies .
Convective-resolving studies have been used for Earth (Rio et al., 2010) and Mars (Colaïtis et al., 2013) atmospheres to improve the parametrization of the convection in the GCMs, the present study could be the starting point of such methodology for a more sophisticated approach for the mixing of heat, momentum, and species in the Venus cloud layer.
A SO 2 Maps