Venus boundary layer dynamics: eolian transport and convective vortex

Few spacecraft have studied the dynamics of Venus' deep atmosphere, which is needed to understand the interactions between the surface and atmosphere. Recent global simulations suggest a strong effect of the diurnal cycle of surface winds on the depth of the planetary boundary layer. We propose to use a turbulent-resolving model to characterize the Venus boundary layer and the impact of surface winds for the first time. Simulations were performed in the low plain and high terrain at the Equator and noon and midnight. A strong diurnal cycle is resolved in the high terrain, with a convective layer reaching 7 km above the local surface and vertical wind of 1.3 m/s. The boundary layer depth in the low plain is consistent with the observed wavelength of the dune fields. At noon, the resolved surface wind field for both locations is strong enough to lift dust particles and engender micro-dunes. Convective vortices are resolved for the first time on Venus.


Introduction
The interaction between the surface and the atmosphere is a major aspect of exchanges of heat and angular momentum, impacting the thermal and wind shear profiles and, therefore, the atmospheric dynamics and the rotation of the solid body itself. On Venus, however, the first 10 km above the surface remain a largely unknown region due to the technological difficulty in probing this region below the cloud layer.
Only a limited number of probes have been able to collect data on Venus. The VeGa-2 probe has successfully measured the only temperature profile in that region (Linkin et al., 1986). At the surface, only Venera 9 and 10 directly measured the wind for respectively 49 min and 90 s (Avduevskii et al., 1977), and several other probes like Venera 13 and 14 measured indirectly the wind speed (Ksanfomaliti et al., 1983). The amplitudes of the measured wind speeds are less than 2 m s´1 below 100 m (Lorenz, 2016), with a higher probability for values below 0.5 m s´1. The height of the planetary boundary layer (PBL) and the diurnal cycle is not known, nor are the effect of the topography. Dunes have been observed in radar measurements with Magellan (Greeley et al., 1992), although the lack of knowledge about the spatial and temporal distribution of winds complicates the interpretation of dust transport.
The Institut Pierre Simon Laplace (ISPL) Venus General Circulation Model (GCM) simulations showed the diurnal cycle of the PBL activity is correlated with the diurnal cycle of surface winds . They observed downward katabatic winds at night and upward anabatic winds during the day along the slopes of high-elevation terrains, resulting in a deeper PBL depth at noon. Yamamoto (2011) performed turbulent-resolving simulations of the Venus PBL with a resolution of 100 m without radiative processes and a variety of near-surface theoretical thermal structures. This experiment yielded a PBL depth below 2.5 km. Morellina and Bellan (2022) studied the turbulent chemical-species mixing at the surface, where high-density-gradient magnitude regions are formed with larger gradients due to the supercritical conditions.
The PBL of the Earth and Mars have been studied extensively in both observations and modelling studies. On Earth, the turbulent flux is in the energy budget of the convective layer. The opposite is true on Mars. For Venus, this budget is poorly known. On Earth, the PBL dynamic is associated with the presence of the water cycle and latent heat release. The understanding of such dynamics is crucial for the clouds and water cycle and the energy balance of the surface (Garratt, 1994). On Mars and Venus, the latent heat is negligible. The depth of the Martian surface convective layer is larger than on Earth, with also a greater temperature diurnal cycle (Hinson et al., 2008). The surface slope winds have a strong impact on the thermal structure of the PBL (Spiga et al., 2011).
In this study, we use Large Eddy Simulation (LES) models developed for the clouds convective regions and to simulate the PBL convective activity at two different locations on the surface. This is done to examine the influence of the topography on the boundary layer convection, and two local times, noon and midnight, to quantify the diurnal cycle.
For the first time, the radiative processes are taken into account in the study of the Venus PBL with prescribed solar and IR heating rates. An additional rate representing the effect of the general circulation, heating/cooling due to large-scale wind advection, is also prescribed.
In Section 2, the model is described. The spatial and local time variabilities of the PBL are discussed in Section 3. The PBL of Venus is compared to the Earth and Mars in Section 4. The impact of PBL characteristics on eolian transport is discussed in Section 5. Our conclusions are summarized in Section 6.

Dynamical core
The LMD LES model 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 over a specified 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 atmospheric variables (winds and potential temperature). The parametrization of the unresolved smallscale eddies is carried out by a subgrid-scale prognostic Turbulent Kinetic Energy closure by Deardorff (1972). This methodology has been used for extensively for Earth convection study (Moeng et al., 2007), and the Martian atmosphere (Spiga et al., 2010), Venus cloud convective layer (Lefèvre et al., 2017(Lefèvre et al., , 2018, and terrestrial exoplanet convection (Lefèvre et al., 2021).

Model Physics
Due to the constant heat capacity of the dynamical core and the architecture of the coupled radiative transfer, as well as computational time, the radiative forcing is handled in the same way that in Lefèvre et al. (2017). The solar and radiative heating rates are extracted from ISPL Venus GCM (Garate-Lopez and Lebonnois, 2018), which uses IR transfer (Lebonnois et al., 2015) based on Eymet et al. (2009), with the latitudinally-varying cloud model of Haus et al. (2014Haus et al. ( , 2015. An additional heating rate is prescribed, representing the large-scale heating from the dynamics. In this study, we focus on the first 10 km and the large-scale heating will come mainly from the anabatic/katabatic slope flows. No surface and sub-surface physics are considered in this study. Lebonnois et al. (2018) showed with GCM modelling that the convective depth in the PBL was impacted by the diurnal cycle of the surface wind, and was maximal in the steepest slope of the equatorial topographic features. Therefore, we choose two locations at the surface with two distinct elevations and slope environments at the Equator. Here, the incoming solar flux is maximized in order to study the activity of the PBL where it is supposed to be the most active. One of the locations is in the low plain, with an elevation of -320 m at 0˝longitude, and the other is in the western part of Ovda Regio with an elevation of 1030 m at 80˝longitude. The point in the plain will be hereinafter referred to as low plain, and the point in Ovda Regio will be referred to as high terrain. The domain of the LES simulations is flat. Due to computational constraints, simulations of an entire Venus day were not possible, and two local times are considered in this study: noon and midnight. The surface is heat flux is set at 90 W m´2 at noon and -1 W m´2 at midnight for the two locations . This flux is constant in time during the simulations over the entire domain, there is no feedback of the PBL turbulence on the sensible flux. For the two locations, the horizontal resolution and timestep are set at 50 m and 0.4 s. However, the size of the surface area varies depending on local time and location, 30ˆ30 km for the high terrain case at noon and 20ˆ20 km for the rest. The vertical resolution and extent also depend on the location and local time, from 10 km above the local surface with a mean resolution of 90 m for the High terrain case at noon to 5 km above the local surface, with a resolution of 60 m. The different horizontal domains size were chosen to allow several connective cells in each horizontal direction, and were determined by trial and error. The different vertical domains size were chosen to allow several kilometers above the convective layer and were based on the Venus IPSL GCM results . To avoid the spurious reflection of gravity waves propagating upward on the top of the model, a Rayleigh damping layer is applied over the last 500 m with a damping coefficient of 0.01 s´1. The heat capacity is set to a constant value over the whole domain of 1181 J K´1, a reference value from the Venus International Reference Atmosphere (Seiff et al., 1985). Fig 1 shows the initial profiles of the temperature, potential temperature, and the different heating rates. The temperature is colder for the high terrain cases, and the diurnal cycle of the temperature is below 3 K for the two location cases. At noon, a neutral layer corresponding to the convective layer is visible below 2 km above the local surface for the low plain case, and below 8 km above the local surface for the high terrain case. At midnight, there is no visible neutral layer, meaning that the convective activity is weak. Regarding the heating rates, the short wave heating is slightly greater for the low plain case, but the thermal cooling is slightly stronger at midnight. However, there is a strong difference for the large-scale heating whereas for the high terrain case it is positive up to 6 km above the local surface, with stronger values at noon. While for the low plain case, the large-scale heating is negative in the first 1.5 km at noon and 2 km at midnight, and then alternating between positive and negative values above. This variability of the large-scale heating reflects the effect of the topography and the diurnal cycle of the surface wind. This variability reflects on the total heating rate, alternating between negative and positive values that will enforce the convective depth. The model is initialized with thermal profiles, winds and radiative rates that reached GCM equilibrium, using hypotheses from the subgrid-scale parametrization that will impact the equilibrium state of the region. With the very poor knowledge of the first 12 km, it is difficult to assess the realism of the initial state used for the present studies, although the surface temperature and winds are consistent with measurements. The time step of the simulations is set to 0.4 s, and the radiative transfer is called every 100 dynamical steps, ensuring that there is no issue with numerical precision (Rafkin and Soto, 2020). The outputs shown in the following sections are obtained after at least 1 Earth day of simulations for the midnight cases, where a steady-state is reached. For the noon cases, after 2 Earth days, there is a small temperature drift due to the set-up of the surface flux and radiative rates, around 10´8 K/s near the surface. A constant forcing corresponding to a single time of day is sometimes used to study self-aggregation in Earth's tropics (Daleu et al., 2015;Wing et al., 2017), as well as tidally-locked rocky exoplanets (Zhang et al., 2017;Sergeev et al., 2020;Lefèvre et al., 2021). When no aggregation is present, the model reaches equilibrium in a couple of days (Wing and Emanuel 2014). The Venus atmosphere can be considered dry. In the present study, there is no water vapor and clouds, and therefore no feedback between moisture and radiation/surface temperature. No aggregation of convection is expected. There is no noticeable change in pressure over time in the domain for all the cases considered. The surface wind's amplitude range does not vary in time. The pseudo-equilibrium is realistic enough to provide qualitative insight into the Venus PBL dynamics.

Simulation settings
The model configuration is summarized in 3 Spatial and temporal variability of the PBL Fig 2 shows snapshots of the vertical and horizontal cross-sections of the vertical wind for the high terrain at noon and midnight, and the low plain at noon. The diurnal cycle for the high terrain location is striking: at noon the elevation of the PBL can reach 7 km above the local surface with vertical wind speed as high as 1.3 m s´1. Whereas at midnight, the convective layer barely reaches 0.5 km with vertical wind speed ă 0.2 m s´1. This diurnal cycle is consistent with GCM simulations . The difference in depth leads to a difference in convective cell diameter. At noon, the typical cell diameter is around 5 km, with some cells reaching 10 km, and about 2 km at midnight. The cellular features are elongated in the ydirection at midnight. The Richardson number is stronger at midnight than at noon, meaning that the shear becomes relatively strong compared to buoyancy. Such features depend on the vertical shear of the horizontal wind, where there is very little data to constrain this result. With similar surface flux, the elevation of the PBL depth is different at noon, 2 km above the local surface for the low plain compared to 7 km for the high terrain. This difference is due to the large-scale forcing: in the high terrain this forcing is positive up to 6 km above the local surface (Fig 1), bringing energy to the environment and leading to a deeper convective layer. Whereas in the low plain, the large-scale forcing is strongly negative in the first 1.5 km thus stabilizing the atmosphere and leading to a smaller PBL depth. In the high terrain at night, the large-scale forcing will also bring more heat, (although less than at noon) to the environment, but the IR cooling is stabilizing the atmosphere. In the low plain at midnight (not shown here), both the large-scale forcing and IR cooling stabilize the atmosphere, resulting in no convective activity. This PBL spatial variability is also consistent with the GCM simulations .
To demonstrate the impact of the large-scale forcing on the convective depth, simulations were carried out without large-scale forcing and wind shear, but with the same thermal profiles. Results are shown in Fig 3. Without the large-scale forcing, the height of the PBL decreases by several kilometers in the high terrain and by several hundred meters in the low plain. At noon, the large-scale circulation will heat the atmosphere. Such an impact on the PBL depth due to the topography is not observed on Earth, but is on Mars (Spiga et al., 2010). The convectionresolving model of Yamamoto (2011) displays a realistic surface heat flux with a PBL depth of Figure 2: Snapshots of the vertical cross-section (left column row) and horizontal cross (right column) at 100 m above the local surface of the vertical wind (m/s) for the high terrain case at noon (top row) and at midnight (middle row) and the low plain case at noon (bottom column).
around 2 km, consistent with the low plain case. With higher surface flux, though unrealistic, the convective layer can go as high as 6 km. In 3D GCMs or 1D models, the turbulence processes are not resolved, therefore a parametrization is needed and a vertical eddy diffusivity coefficient, or K zz , is often used to represent convection. In LES models, like the one presented here, the resolution is small enough to resolve the larger eddies, and no vertical eddy diffusivity is therefore needed.
K zz has been estimated in the Venus atmosphere modelling, but is still not well-constrained. The 1D models of Matsuda and Matsuno (1978) and Takagi et al. (2010) have estimated the vertical eddy diffusion between 10 2 and 10 3 m 2 s´1. In the IPSL Venus GCM, the subgrid parametrization exhibits values as high as 1 m 2 s´1 in the high terrain . Yamamoto (2011) convection-resolving model displays values between 1 and 10 4 m 2 s´1.
Using Prandtl mixing-length theory (Lindzen, 1971) the vertical eddy diffusion can be estimated from the resolved convective plume and, as expected, it is higher with a stronger convective activity. It reaches values as high as 10 3 m 2 s´1 in the high terrain at noon, consistent with the estimation from Yamamoto (2011) with similar PBL depth. However, the values are several orders of magnitude higher than the IPSL Venus GCM. The parametrization of the GCM, based on Mellor and Yamada (1982), does not display realistic values for vertical eddy diffusion compared to the convective-resolving model of Yamamoto (2011) and the present study.

Comparison with Earth and Mars
The depth of the convective layer in the low plain at noon is comparable to the Earth's shallow convection depth (Garratt, 1994) but with much weaker wind, 0.2 m s´1 for Venus, and several meters per second for Earth (Park et al., 2018). The depth of the convective layer in the high terrain is somewhat comparable to the Martin PBL depth, reaching 10 km at maximum (Hinson et al., 2008). This has a much lower vertical wind amplitude, around 1 m s´1 for Venus, and 10 m s´1 for Mars (Spiga et al., 2010). The difference in vertical wind is also visible in the convective velocity scale, defined as with z h the height of the PBL. For Venus, W˚is around 0.4 m s´1, around 2 m s´1 for Earth (Stull, 1988), and between 4 and 6.5 m s´1 for Mars (Spiga et al., 2010). Further comparison between the three atmospheres is shown in Fig 4 with the vertical eddy heat flux and vertical velocity variance in the PBL in dimensionless defined respectively as The heat flux on Earth is maximum at the surface, whereas for Venus and Mars fluxes are at the surface, and decrease with height. The turbulence acts to heat the atmosphere. On Mars and Venus, the turbulent flux first increases with height, acting to cool the atmosphere, and then decreases with height, as on Earth, to heat the atmosphere. The low vertical wind amplitude in the PBL shows that with a denser atmosphere, around 50 times higher than the Earth, and 3000 times that of Mars atmospheric density, the heat transport by the atmosphere is more efficient than on other terrestrial planets. For Venus, as on Earth, the energy budget of the convective layer is dominated by the convective flux, with negligible solar and IR heating rates. In the Martian atmosphere, the radiative heating rates are dominant, with surface temperatures that do not depend on altitude, and a pressure effect that is obtained for similar solar heating and IR cooling rates (Spiga et al., 2010). The empirical similarities in this dimensionless approach between the three planets provide a framework for boundary layer parametrization. Figure 4: Vertical profiles of the vertical eddy heat flux (left) and vertical velocity variance (right) for Venus, Earth and Mars PBL. Heat flux, variance and height are normalized respectively by the maximum vertical eddy heat flux, the convective velocity scale W˚and the boundary layer depth. For Venus, the high terrain case at noon is presented, Earth data is from Stull (1988), and Mars data is from Spiga et al. (2010).

Eolian Transport
Only two prominent dune fields have been identified on Venus with the Magellan radar (Greeley et al., 1992), Algaonice at 25˝S, 340˝E cover some 1300 km 2 with bright dunes of 0.5-5 km in length with a wavelength around 0.5 km, and Fortuna-Meshkenet at 67˝N, 91˝E with transverse dunes of 0.5-10 km long, 0.2-0.5 km wide and spaced by an average of 0.5 km. Several sites with micro-dunes or small-scale wind streaks have also been proposed (Weitz et al., 1994;Bondarenko et al., 2006) but the resolution of the Magellan radar is not high enough to confirm them.
The wavelength of the observed dune fields is about 500 m. The height of the PBL in Algaonice dune field would be slightly lower than the one at the Equator (Fig 2-bottom line) with an averaged depth over a Venus day around 500 m, about the same order of magnitude than the wavelength of the dunes. However, in the Fortuna-Meshkenet the averaged PBL depth would be much lower due to the solar flux more than ten times lower, although the proximity of Ishtar Terra could engender wind turbulence. On Earth, the wavelength of giant dunes is controlled by the depth of the PBL (Andreotti et al., 2009). On Venus, From the two Magellan dune fields observations, the relationship between the wavelength of a giant dune and the depth of the PBL seems to be more complicated, and more observations are needed.
The wind profile is generally approximated by the Prandlt-von Kármán equation defining the friction velocity u ‹ as: with the Von Kármán constant, k, equals to 0.4, the surface roughness, z 0 , set to 1 cm. u z is the horizontal wind at the altitude z. Theoretical calculations and laboratory experiments estimated the threshold friction velocity for which the dust is lifted in Venus surface condition to be minimum around 2.5 10´2 m s´1 (Iversen et al., 1976), about 10 times lower than for Earth, and 100 times lower than for Mars due to the very dense atmosphere, depending on dust radius. Such a value depends also on the density of the dust, data that is not well-known yet. Laboratory experiments tested the dust transport in Venus surface conditions and showed that for surface velocities between 0.63 and 1.5 m s´1 there was the formation of micro-dunes, with a wavelength between 8 and 27 cm (Greeley et al., 1984). Below 0.63 m s´1, the amplitude of the wind is too small to transport dust, and above 1.5 m s´1 no dunes are formed because dust grains can be transported across the lee of a wave onto another and therefore blurring the separation between dunes.
For both the high terrain and low plain, both at noon, the friction velocity of the vast majority of the spatial area is above the saltation threshold, meaning that the turbulence is able to transport dust particles, as shown in Fig 5. On the other hand, at night for both locations, the friction velocity is below the saltation threshold over the entire area, and no dust is transported. At night, the horizontal wind amplitude distribution at 10 m above the local surface for both locations is fully consistent with in-situ measurements (Lorenz, 2016), with wind speed ă 0.5 m s´1. Whereas at noon, the majority of the wind distribution is ą 0.5 m s´1, with values of 1 m s´1 for the high terrain. These values are on the outer edge of the observed values, without taking site-specific factors like slopes or local times into account. With a deeper convective layer in the high terrain, the distribution is broadened compared to the low plain location. In the low plain, the horizontal wind amplitude is almost always in the amplitude range where micro-dunes are forming, whereas in the high terrain a significant part of the horizontal wind is too weak to transport wind to form micro-dunes. The formation of micro-dunes is therefore more probable in the low plain at noon, however the spatial distribution of particle reservoir is not known. Fig 6 shows maps of the surface temperature and surface pressure in the high terrain at noon. The temperature anomaly induced by the convection is very small, around 0.05 K for the strongest convective activity. However, in the pressure field, there is a significant drop, characteristic of convective vortices at the convergence of convective cells .
At these convergences, a pressure deficit is formed from rising plumes, where the updrafts cause a vertical stretching of horizontal wind, leading to a spiralling flow of rising air (Sinclair, Figure 5: Snapshots of the friction velocity u ‹ (top row) and horizontal wind amplitude distribution (bottom row) in the high terrain (left column) and low plain (right column) at 10 m above the local surface, both at noon. 1973; Rennó et al., 1998;Toigo et al., 2003;Ito et al., 2013). The vorticity inside the vortices, visible with the vertical component of the rotation vector ζ = (Bv/Bx) -(Bu/By) with u and v the two horizontal wind fields. There are vortices with both positive and negative vortices, meaning, respectively, cyclonic and anticyclonic advection. These patterns are also visible in the horizontal divergence D = (Bu/Bx) + (Bv/By), with convergence (positive divergence values) for the cyclonic advection and divergence (negative divergence values) for anticyclonic advection. The vortices are also regions where the shear, Bu/By + Bv/Bx, and stretching, Bu/Bx -Bv/By, deformations are strong. Fig 7 shows horizontal cross-sections of the surface pressure, vorticity, divergence, shear and stretching at y = 25.2 km (black line in Fig 6) where there is the presence of a strong vortex. The main vortex is visible at x = 10 km, with the presence of two smaller vortices at x = 12.2 km and x = 12.9 km with lower pressure drop. The vortices are strongly visible in the vorticity field, with values up to 2ˆ10´2 s´1, at least an order of magnitude lower than Earth vortices (Spiga et al., 2016).The vortices are barely discernible in the divergence fields, but generate a lot of shear and stretching. There are generations of several strong shear and stretching increase or drops, without noticeable pressure decreases at x = 14.7 km, 22 km or 23.9 km, for example, due to the presence of thin convective updrafts between convective cells. The main vortex extends up to 5 km above the local surface, much higher than on Earth, where it is capped by the height of the PBL (around 2 km). The two smaller vortices only extend up to 1 km above the local surface. The lifetime of such vortices can be several hours. No such structures are resolved in the low plain at noon with a 50 m resolution. The convective activity is either too weak to form vortices, or the resolution of the simulations is too coarse to resolve thin vortices. At night for both locations, no vortices are resolved because the convective activity is too weak to create a vortex. The pressure drops are around 100 Pa for the largest vortices, comparable to the highest pressure drops of dust devils on Earth (Lorenz and Lanagan, 2014). However, the diameters are different, almost 1 km on Venus against a few hundreds of meters, maximum, for Earth dust-devils (Lorenz, 2011). Convective vortices are also measured in the Martian atmosphere (Baker et al., 2021), with pressure drops being inferior to a few pascals Spiga et al., 2021). The 50 m horizontal resolution used for these simulations is sufficient to qualitatively resolve the characteristic structure of convective vortices, but a horizontal resolution one order of magnitude lower would be necessary to resolve a larger part of the spectrum of kinetic energy (Spiga et al., 2016), as well as smaller vortices. Lorenz (2021) used an equation to estimate the maximum vertical wind amplitude in a convective vortex for the different terrestrial planets and Titan. For Venus, this maximum vertical wind amplitude was estimated at 2 m s´1, stronger than the one in the present study around 1.3 m s´1. However, the value of the heat capacity and surface sensitive heat flux chosen in Lorenz (2021) were not representative of the Venus surface at the Equator. With a heat capacity of 1181 J K´1 and a surface sensitive heat flux of 90 W m´2, the estimation with the equation from Lorenz (2021) is 1.35 m s´1, very close to the value obtained in the present model in the high terrain. The location of the vortices is where the friction velocity is the highest and superior to the saltation threshold. Therefore, if there is dust material to be lifted, it will generate dust devils.

Conclusion
This study presents turbulent-resolving simulations of the Venus PBL, with for the first a realistic thermal profile and with the impact of the large-scale circulation. Simulations were performed at two locations on the Equator, in the low plain and high terrain, and at two local times, noon and midnight. The convective activity is stronger at noon due to the incoming stellar radiation. It is also stronger in the high terrain due to the large-scale heating from the slope winds. The vertical eddy diffusion estimated from turbulent-resolving simulations is not consistent with the GCM, with values several orders of magnitude higher, consistent with previous turbulent-resolving modelling. The PBL of Venus was compared to its equivalent in the atmosphere of the Earth and Mars. On Venus and Mars, the turbulent flux first increases with height, acting to cool the atmosphere, whereas on Earth the turbulence acts to heat the atmosphere. The energy budget of the convective layer is dominated by the convective flux, with negligible solar and IR heating rates for Venus and the Mars. In the Martian atmosphere, the radiative heating rates are dominant. At noon, the turbulent activity is strong enough for the friction velocity to be above the saltation threshold needed for eolian transport. However, there are a lot of unknowns about the composition, density, radius, and reservoir of such particles that could drastically affect their transport. These simulations were carried out at noon on the Equator, where the convective activity is more intense. More simulations need to be performed at higher latitude and other local times. In the high terrain at noon, the strong activity in the PBL convective vortices are resolved for the first time. Pressure drops as strong as 100 Pa are observed, comparable to the highest pressure drops of dust devils on Earth, but with a larger diameter, almost 1 km on Venus versus a few hundreds km for the Earth. The vertical velocity inside such vortices can exceed 1 m s´1. With the presence of dust, these vortices could lead to dust devils.
One of the main improvements of the model setup would be to have an interactive surface heat flux and radiative transfer like in the Venus cloud layer (Lefèvre et al., 2018), this would allow feedback between the surface and the atmosphere.
In the IPSL Venus GCM, a nitrogen vertical gradient in the first 7 km from the surface has been tested (Lebonnois and Schubert, 2017;Lebonnois et al., 2018) and has a substantial impact on the stability of the atmosphere in this region and leads to weaker PBL depth. Such gradient should be tested with a similar setup in future studies.