Analyzing the Instabilities in the Venus Atmosphere Using Bred Vectors

Barotropic, baroclinic, and Rossby‐Kelvin instabilities exist in the Venus atmosphere, as revealed by previous studies using numerical models. This study aims to deepen our understanding of the instabilities in the Venus atmosphere using the breeding of growing modes (BGM) method for the first time. We conducted a 1‐year model simulation as the control run. First, we conducted identical twin experiments in which random perturbations were added to the control run at the initial time, and they grew freely. Perturbations in the upper layer (UL, 70–100 km altitudes) grow faster at the beginning but saturate earlier at a low value compared to those in the lower layer (LL, 40–70 km altitudes). Next, we implemented the BGM method and conducted multiple breeding cycle experiments with different rescaling amplitudes and intervals, the two parameters of the BGM method. The bred vectors (BVs) from these experiments identified instabilities in various regions. With large rescaling amplitudes, the structure and amplitude of the BVs in the LL closely resembled the deviations in the control run, indicating the growth of BVs due to barotropic, baroclinic, and Rossby‐Kelvin instabilities. Composite mean analysis shows larger BV amplitudes in the morning hemisphere at the 60–70 km altitudes in the mid‐latitudes, indicating enhanced baroclinic instability by the thermal tide. Finally, we estimated the intrinsic predictability related to baroclinic instability for Venus is >1 month, which would be longer than that for the Earth of ∼2 weeks.


Introduction
Numerical simulation plays an important role in understanding the Venus atmosphere.The Atmospheric General Circulation Model (AGCM) for the Earth Simulator (AFES) (Ohfuchi et al., 2004) was adapted for the Venus atmosphere, resulting in the AFES-Venus model (Sugimoto et al., 2014a(Sugimoto et al., , 2014b)).The AFES-Venus model successfully simulated essential features of the Venus atmosphere, including the superrotation and the thermal tide (Sugimoto et al., 2014b;Suzuki et al., 2022;Takagi et al., 2018).Furthermore, simulations using this model have revealed various instabilities in different regions (Ando et al., 2017;Sugimoto et al., 2014aSugimoto et al., , 2014b;;Takagi et al., 2022).In the high latitudes, simulations showed that the zonal wave number one component of the axiasymmetric temperature disturbance is predominant at the 50-75 km altitudes (Ando et al., 2017).The vertical structure of these disturbances is consistent with radio occultation measurement, and they can be interpreted as neutral barotropic Rossby waves associated with barotropic instability in the polar region.In the mid-latitudes, simulations revealed baroclinic instability at the cloud level, where the superrotation with vertical shear is dynamically balanced with temperature (Sugimoto et al., 2014a).In the low latitudes, simulations successfully reproduced Kelvin waves (Takagi et al., 2022).The Kelvin waves are coupled with the Rossby waves in the mid and high latitudes at critical latitudes, and they are excited by Rossby-Kelvin instability.
The breeding of growing modes (BGM) method (Toth & Kalnay, 1993) can elucidate the physical origin of the instabilities and help improve ensemble forecasting and data assimilation.The BGM method obtains growing perturbations from "Breeding Cycle" (BC) experiments.In chaotic systems, two states that are initially similar diverge over time, and the difference between them is known as the "bred vector" (BV) in the BGM method.BV is related to the physical instabilities in the model (Greybush et al., 2013;Toth & Kalnay, 1993).While the growing modes in BV amplify due to instabilities, the decaying modes lose their amplitudes (Toth & Kalnay, 1997).Among the growing modes, for the Earth's atmosphere, the fast-growing modes saturate earlier at lower amplitudes (Figure 1a red solid line).In comparison, the slow-growing modes saturate later at higher amplitudes (Figure 1a black solid line) (Toth & Kalnay, 1997).
The breeding cycle has two parameters: the rescaling amplitude (R amp ) and the rescaling interval (R int ).Smaller R amp and R int values usually emphasize the fast-growing modes, whereas larger R amp and R int highlight the slowgrowing modes.The role of these two parameters can be demonstrated using Figure 1a.When R amp (red dash line) is smaller than the saturation amplitude of the fast-growing modes, the fast-growing modes dominate the total amplitude at the end by growing faster than the slow-growing modes.Conversely, if R amp (black dash line) is set higher than the saturation amplitude of the fast-growing modes, the fast-growing modes stop growing upon reaching saturation.In contrast, the slow-growing modes can continue their growth and eventually dominate the total amplitude.For the Earth's atmosphere, small R amps emphasize convective instability, while large R amps emphasize extratropical baroclinic instability (Toth & Kalnay, 1993, 1997).In addition, Peña and Kalnay (2004) constructed an idealized model based on the Lorenz 63 model and incorporated different instabilities with different timescales.They showed that the breeding cycle experiments successfully extracted different growing modes by setting different R amp and R int parameters.
The BGM method can elucidate the physical origin of the instabilities.The energy equations for the BV were developed to diagnose the baroclinic and barotropic energy conversions for the tropical instability waves in the tropical Pacific Ocean (Hoffman et al., 2009).These equations were also used to explain the physical origins of the instabilities in the Mars atmosphere (Greybush et al., 2013).The dashed lines represent the rescaling amplitudes (R amps ).This figure is adapted from Figure 7 in Toth and Kalnay (1993).(b) Steps involved in conducting breeding cycle experiments to obtain bred vectors (BVs).
The BGM method is also valuable for improving ensemble forecasts and data assimilation techniques.To capture the true evolution of the atmosphere by ensemble forecasts, it is necessary to perturb the analysis along the growing errors present in the analysis (Toth & Kalnay, 1993).Achieving this involves running multiple breeding cycles, each starting from distinct initial conditions obtained by adding different perturbations to the initial condition of a control run.Due to nonlinear interaction, stochastic forcing, and regional features, BVs from different breeding cycles span the subspace of the growing perturbations (Toth & Kalnay, 1997).These BVs are added to the analysis to initialize ensemble forecasts.They can also be used to calculate the bred vector dimension (BVD), which measures the effective local dimensionality of a chaotic system (Patil et al., 2001).
To the best of the authors' knowledge, applying the BGM method to investigate instabilities in the Venus atmosphere has not been explored thus far.In this study, we use this method to gain new knowledge of the instabilities of the Venus atmosphere.Our investigation focuses on several aspects: (a) analyzing the amplitudes and growth rates of perturbations influenced by different instabilities, (b) comparing the BVs with deviations from the zonal mean in the control run, and (c) examining the impact of the thermal tide on the instabilities and BVs.This paper is organized as follows: Section 2 introduces the data and method.Section 3 presents the results obtained from the breeding cycle experiments.Section 4 addresses relevant issues.Finally, Section 5 provides a conclusion.

Model Setting
Our setting of the AFES-Venus model is the same as that of the data assimilation (DA) system for the Venus atmosphere (Fujisawa et al., 2022;Sugimoto et al., 2017) because we intend to understand further the dynamics of the model used in the DA system.The setting also resembles the setting in Sugimoto et al. (2014b).Although Venus rotates westward, we set its rotation to eastward to compare the dynamics of the Earth's atmosphere.The primitive equations in sigma coordinates are used in the model, but the surface is assumed to be flat.Observational studies (Schubert et al., 1980) and model simulations (Imamura et al., 2014) suggest convective activities at the cloud level.However, as our model is hydrostatic, it cannot treat the convection directly.To represent the convection implicitly, we use the convection adjustment scheme (Sugimoto et al., 2014a).In this scheme, the temperature lapse rate is restored to neutral when an atmospheric layer becomes statically unstable.The triangular truncation number of the model is set to 42 (T42), resulting in a resolution of 128 (zonally) × 64 (meridionally) × 60 (vertically).The model time step is 150 s.The model reaches an altitude of 120 km with a vertical grid spacing of 2 km.The model does not incorporate moist processes and the planetary boundary layer.Horizontal and vertical eddy diffusions are included.The horizontal eddy viscosity is represented by second-order hyperviscosity, with a damping time of approximately 0.1 Earth days for the largest wavenumber components.The vertical eddy diffusion coefficient is set at 0.15 m 2 s 1 .Rayleigh friction is implemented at the lowest level to mimic surface friction.Additionally, a sponge layer is introduced above 80 km altitude solely for eddy components.Solar heating is applied based on the vertical and horizontal distributions outlined in Tomasko et al. (1980).Newtonian cooling is used for the radiative forcing in the infrared region, with coefficients based on Crisp (1989).The temperature is relaxed to a reference temperature that is horizontally uniform based on observations.The model output includes zonal wind (ms 1 ), meridional wind (ms 1 ), temperature (K), surface pressure (Pa), etc.The reader can refer to Sugimoto et al. (2014b) for more detailed information on the model settings.Previous studies about the instabilities applied slightly different model settings because they wanted to study some particular phenomenon.For instance, Ando et al. (2017) considered solar heating below 90 km, while in our study, it is below 80 km.In Takagi et al. (2022), the static stability and solar heating profiles differ slightly from this study.With these differences in mind, we will compare our findings to theirs.
To explore the impact of the thermal tide on instabilities, we designed two groups of experiments.Given that the diurnal component of solar heating plays a key role in generating the thermal tide (Sugimoto et al., 2017), we conducted experiments that included (Qt) and excluded (Qz) the diurnal component.For each group, we performed multiple breeding cycles using various combinations of R amps and R ints .In this paper, we primarily focus on presenting the results from the Qt group.However, we also compare the results from the Qz group with those of the Qt group when significant differences arise.

Definition of the BV Norm
To compare the amplitude of different modes in the perturbation and to appropriately rescale BV in the breeding cycle, it is necessary to define a norm for the perturbation.Generally, BVs are not sensitive to the choice of the norm (Corazza et al., 2003;Kalnay, 2002;Toth & Kalnay, 1997), and different norms have been utilized in various studies, such as stream function, temperature squared, and others.However, Pazó et al. (2011) argued that choosing a norm can influence the BVs, as demonstrated in their study.We propose that a reasonable norm should allow for a fair comparison of amplitudes between different modes.Toth andKalnay (1993, 1997) demonstrated that the low-energy convective mode and the high-energy baroclinic mode can be emphasized by using different R amps and R ints .Given that the evolution of BVs involves energy conversion, we consider the energy norm as a suitable candidate for this purpose.
The total energy norm has been used in several studies (e.g., Magnusson et al., 2008, their Equation 1).It includes the contribution from the surface pressure squared term.Since our objective is to compare the amplitudes of modes at various latitudes and heights, if we use the same total energy norm, we encounter the challenge of assigning the contribution of the surface pressure squared term to the amplitudes of the perturbations at different heights.Consequently, our definition of the energy norm excludes the surface pressure squared term, allowing for a more consistent assessment of the amplitudes of the perturbations at different latitudes and heights.In comparison, some studies (e.g., Thanh et al., 2016;Zhang et al., 2003) also ignored the surface pressure squared term from their total energy norm, although their total energy norm is the total energy per unit mass.Our definition of the energy norm can be expressed as follows: where δu (ms 1 ), δv (ms 1 ), and δT (K) represent the zonal wind, meridional wind, and temperature perturbations, respectively.dL (m) corresponds to the longitudinal distance of a model grid at the equator, while dy (m) represents the latitudinal distance of a model grid.g denotes the gravitational force constant, with a value of 8.87 m s 2 .c p represents the specific heat at constant pressure, with a value of 1,000 Jkg 1 K 1 .T r denotes the reference temperature, with a value of 737.15 K, and p s represents the surface pressure (Pa) in the control run.dσ refers to the difference between two sigma levels, and θ corresponds to the latitude.Σ x Σ y Σ z indicates horizontal and vertical summation.The first two terms in the expression correspond to the kinetic energy of the perturbation, while the third term is part of the available potential energy of the perturbation (Keller et al., 2008).The unit of E is joule (J).Since the term dLdy 2g is constant for every grid point and its value is large, we divide E by a constant number 5*10 17 to facilitate calculation and visualization.At each grid point, we calculate where its unit is J, with the normalization factor 5*10 17 .Various regions of summations of E grid are used for different purposes.To rescale the BV in the breeding cycle, the total amplitude of the perturbations is computed by summing E grid from 10 to 110 km altitudes and horizontally across the entire globe.This range is chosen as we primarily focus on the instabilities at these levels.In addition, to compare the amplitudes of the perturbations across different regions, we computed the summation of E grid in different regions.For simplicity, in the subsequent discussion, we will denote the amplitude using numerical values without including the unit and the normalization factor.The actual amplitude value will be the value multiplied by 5*10 17 .For example, if the amplitude value is 1,000, the amplitude equals 1,000 * 5* 10 17 J.If the breeding cycle is conducted using R amp = 1,000 and R int = 6 hr, we name this experiment R amp 1000-R int 06.
Additionally, we explored alternative norms, such as the temperature-squared norm and the total energy norm that includes the surface pressure squared term.Preliminary analysis indicates that the results obtained using these norms are qualitatively similar to the energy norm defined in Equation 1.However, conducting a comprehensive analysis of these alternative norms and exploring additional norm choices are beyond the scope of this study.Therefore, in this study, we solely use the norm defined in Equation 1 to perform the breeding cycle and quantify the amplitudes of the perturbations in different regions.

Breeding Cycles
To obtain the BVs, multiple breeding cycles were performed using different values of R amps and R ints .Section 3 will demonstrate how to set these two parameters to emphasize different instabilities.The following steps outline the procedure for conducting the breeding cycle.The step numbers below correspond to the same numbers in Figure 1b.
1. Conducting a control run.In this study, we use the YYYY/MM/DD/HH format to represent time, where, for instance, 0005/01/01/01 indicates the fifth year, first month, first day, and first hr.Also, unless otherwise specified, all time information in this study is based on Earth time.We first initialized the model from an idealized eastward super-rotating flow.The zonal mean zonal wind increases linearly from the surface to 70 km altitude and is constant above.It reaches a maximum value of 100 ms 1 at 70 km at the equator.The meridional distribution of the zonal wind follows solid-body rotation.We assumed a generalized gradient wind balance: the temperature field is balanced with the zonally uniform zonal flow (Sugimoto et al., 2014a).Starting from this initial state, a simulation was conducted for 5 years.The model atmosphere reaches a quasi-equilibrium state within the first year, and the temperature and wind fields remain stable for over 10 years.The model states in the fifth year (from time 0005/01/01/00 to 0006/01/01/00) were selected as the control run for the breeding cycles.2. To obtain the initial condition for the perturbed run, a random perturbation is added to the control run at the initial time (t0: 0005/01/01/00 in this case).The random perturbation is obtained through the following process: For each variable at each grid point, a scalar value is drawn randomly from a Gaussian distribution.
The mean of the Gaussian distribution is set to 0, while the standard deviation is calculated from the time series of a 4-year simulation starting from 0005/01/01/00. 3. The perturbation at the initial time is rescaled to have an amplitude of R amp .The method for rescaling is similar to step ( 6), but regardless of whether the amplitude of the perturbation is larger or smaller than the R amp .After the perturbation is rescaled, it is added to the control run to obtain a new perturbed state (p′(t0)).4. The perturbed run, initiated from the perturbed state (p′(t0)), is evolved for a certain period (R int ) until the next rescaling time (p(t1)).During this period, certain growing modes within the perturbation will exhibit faster growth than others, resulting in their amplitudes gaining a larger proportion of the total perturbation amplitude.5. To obtain the perturbation at the rescaling time, the control run is subtracted from the perturbed run (p(t1) c( t1)).6.The perturbation is rescaled.After conducting multiple cycles, the perturbation is denoted as BV.If the norm of the BV (‖BV‖ 2 ) is smaller than the R amp , we do not rescale the BV.In case when ‖BV‖ 2 is larger than the R amp , BV is rescaled such that its amplitude becomes equal to the R amp .Firstly, the rescaling factor (α) is computed as α = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ Ramp/ ‖BV‖ 2 √ .Secondly, BV is rescaled by multiplying it by α to obtain a new BV, denoted as BV′ (BV′ = BV* α).All components of BV, including the zonal wind, meridional wind, temperature, and surface pressure, are rescaled by the same factor.After rescaling, the direction of the new BV remains the same as the original BV, while its amplitude is changed to the The new perturbed state at the rescaling time is obtained by adding the rescaled BV to the state of the control run (p′(t1) = c(t1) + BV′).
Steps ( 4) to ( 7) are iterated over a certain period.Eventually, the pattern of the BV will stabilize, with certain modes exhibiting larger amplitudes than others.The breeding cycle experiment was conducted for 1 year in this study.Various combinations of R amp and R int were tested.
The growth rate of the BV characterizes the instability.To calculate the daily growth rate of BV, we followed the method outlined in Greybush et al. (2013).Let's denote the amplitude of BV at the current time (t) as N t .It is important to note that if BV has been rescaled at the current time, N t is equal to the R amp .Assuming the norm of BV at the next time step (t + δt) before rescaling is N t+δt , we can compute the growth rate within this interval as G t→t+δt= ln( N t+δt N t ) /δt, where "ln" represents the natural logarithm.In our simulation, data was saved every 6 hr.We calculated the 6-hourly growth rate G 0→6h for the interval from 0 to 6 hr.Similarly, G 6→12h , G 12→18h , and G 18→24h were calculated.The daily growth rate, denoted as G 0→24h , was obtained by summing these individual growth rates ) .In addition, the daily growth rates were averaged to obtain the average growth rate within a specific period we are interested in.By varying the combinations of R amps and R ints , we obtained the daily growth rate and the average growth rate within one Venus year.

Basic State and Zonal Deviations in the Control Run
We present the latitude-height cross sections of the temporal and zonal average of the zonal wind (u), meridional wind (v), and temperature (T) in the control run (Figures 2a-2c contour).The averaging period spans from 0005/ 02/01/00 to 0005/09/15/18, equivalent to one Venus year.The root mean (temporal and zonal average) square of the deviations (u′, v′, T′) from their zonal means are superimposed on the same figures (Figures 2a-2c shading).The wave amplitudes at higher altitudes are generally larger because density decreases dramatically with height.To visualize the waves at lower atmospheres, we multiply the deviations by the square root of the sigma level difference ( ̅̅̅̅̅ dσ

√
).The figures reveal that the wave amplitudes are symmetric about the equator.According to previous model studies, 4-day and 5-day waves at the Venus cloud top are symmetric about the equator (Takagi et al., 2022), while the 7-day wave is antisymmetric about the equator (Takagi et al., 2023).In our study, we did not decompose the waves with different wave numbers.The evolution of the waves is generally asymmetric about the equator.However, we found that the conclusion of this study is also true for the southern hemisphere.Therefore, we primarily show our analyses for the northern hemisphere.To facilitate the examination of waves in different latitudes, we divide the latitudes in the northern hemisphere into three regions: low (0-30°N), mid (30-60°N), and high latitudes (60-90°N) (Figure 2d).In addition, to understand the vertical structure of the waves, the longitude-height cross sections of u′, v′, T′ at 80°N, 45°N, and 0°at time 0005/06/03/18 are depicted in Figure 3.We selected this particular time because the waves at that moment exhibit typical structures related to the instabilities.The reader can observe the time evolution of these waves in the animations in Movies S1-S3.The barotropic, baroclinic, and Rossby-Kelvin instabilities have been investigated by Ando et al. ( 2017), Sugimoto In LL, it is divided into 6 regions with a meridional interval of 30°.Labels A to F correspond to the same labels in Figure 4c.et al. (2014aet al. ( , 2014bet al. ( ), and Takagi et al. (2022)), respectively.Notice that their studies separated fast and slow waves (e.g., thermal tides) and mainly focused on the fast wave.Despite such differences, we still compare their studies to ours (Tables 1-3).
In the high latitudes 60-90°N, specifically in the latitudes larger than 75°N, the zonal wind exhibits weak vertical shear but strong meridional shear around the 50-70 km altitudes (Figure 2a, contour).Correspondingly, the meridional temperature gradient is small in this area (Figure 2c, contour).Such atmospheric conditions create favorable conditions for barotropic instability, which can lead to the generation of Rossby waves, as suggested by Ando et al. (2017).The large magnitude of the deviations at the 40-70 km altitudes in >75°N region indicates the presence of these waves (Figures 2a-2c, shading).Since we did not separate the fast and slow waves like those in Ando et al. (2017), the observed waves in this study include both barotropic Rossby waves and the thermal tide.To understand the vertical structure of these waves, we analyze the longitude-height cross section at 80°N (Figures 3a,

̅̅̅̅̅ ̅ dσ √
).The meridional wind deviations in (d) are superimposed on temperature deviations and shown as contours in (g).(b, e, h) are similar to (a, d, g) but at 45°N.(c, f, shading in i) Are similar to (a, d, shading in g), but at 0°N.However, the zonal wind deviations in (c) are superimposed on temperature deviations (shading in (i)) and shown as contours in (i).
3d, and 3g).During the time evolution (Movies S1a, S1c, and S1e), the zonal wave number one component of the deviations is predominant at the 40-70 km altitudes.In comparison, the vertical range of the deviations is around 50-75 km in Ando et al. ( 2017) (Table 1).The difference could be due to the model setting and the vertical weighting.At the same altitudes, u′ (Figure 3a) and v′ (Figure 3d) have constant phases vertically, and they are zonally out of phase by π/2.T′ has a constant phase vertically at 40-65 km, but its phase at 65-70 km is different by π (Figure S3g in Supporting Information S1).Below the 65 km altitude, v′ and T′ are out of phase by π/2 in the zonal direction, while u′ and T′ are zonally in phase.The latitude-longitude cross sections of the deviations (Figure S1 in Supporting Information S1) at the 47, 55, and 70 km altitudes reveal the dominance of the polar vortex between 60°N and 90°N.At the 57 km altitude, the temperature deviation is within ±6 K in Ando et al. ( 2017) and our study (Table 1).Although our study incorporates the thermal tide, the fast waves remain dominant below 70 km.This is evident in the animations where the waves propagate eastward below 70 km (Movies S1a, S1c, and S1e).Estimated from the animation, the period of these waves is around 7-8 days.It is slightly longer than the numbers in Ando et al. ( 2017), possibly because we did not separate the fast and slow waves.In general, the characteristics of the waves below 70 km in the high latitudes observed in our study are similar to those documented by Ando et al. (2017).They proposed that these waves are the neutral barotropic Rossby waves maintained by barotropic instability at 60-68 km levels in the polar region.Conversely, above 70 km, the thermal tide exerts a strong influence and dominates, resulting in a gradual westward movement of the waves, as demonstrated in the animations.
In the mid-latitudes (30-60°N) at the altitudes of 50-70 km, there is a strong vertical shear of the zonal wind, as depicted in Figure 2a contour.Concurrently, the temperature gradient between the equator and the pole exceeds 25 K at these vertical levels, indicating a pronounced meridional temperature gradient (Figure 2c).Furthermore, the atmospheric layer in this region exhibits weak stratification (not shown).These conditions create a favorable environment for baroclinic instability, which can generate waves, as suggested by Sugimoto et al. (2014aSugimoto et al. ( , 2014b)).While high amplitudes of u′ and v′ can be observed in this region, they are relatively smaller compared to the polar region (Figures 2a and 2b shading).The presence of the thermal tide is notable around the 65-70 km altitudes, as indicated by the large amplitudes of u′, v′, and T′ at this level.Figures 3b, 3e, and 3h demonstrate that zonal wave number one component of u′, v′, and T′ dominate around the 50-70 km altitudes at 45°N.v′ and T′ are slightly out  of phase in the zonal direction, with a tilting structure from lower-east to upper-west.In Sugimoto et al. (2014b), similar structures can be seen at around 50-70 km (Table 2).The temperature deviation at the 60 km altitude ranges from 5 to 5 K in Sugimoto et al. (2014b) and in our study.Generally, the structure, amplitude, and development of the waves in our simulation resemble the findings of Sugimoto et al. (2014b), which propose that the baroclinic Rossby waves in the mid-latitudes are generated by the baroclinic instability.Eastward-moving waves dominate below 65 km, as observed in the animation (Movies S2a, S2c, and S2e), indicating the prevalence of Rossby waves.Conversely, the thermal tide dominates at the 65-70 km altitudes, characterized by westward slow-moving waves in the animation at these levels.
In the low latitudes (0-30°N) at the altitudes of 52-60 km (Figures 2a and 2b  3).Therefore, these waves are associated with the Rossby-Kelvin instability.Lastly, above 60 km, the deviation fields are dominated by the thermal tide, indicated by the westward slow-moving waves in the animation (Movies S3a, S3c, and S3e).

Identical Twin Experiments
To determine the appropriate R amp and R int in the breeding cycle, we analyze how the perturbations grow freely in different regions and assess their saturation amplitudes.This is usually fulfilled by conducting identical twin experiments (Judt, 2018).The following steps indicate how to conduct these experiments.
1. 15 random perturbations were obtained following the method described in Section 2.3 step (2).
2. The perturbations were rescaled to have an amplitude of 1 (J, normalization factor 5*10 17 ), followed the method of step (6) in Section 2.3.3. 15 perturbed states were obtained by adding these rescaled perturbations to the state of the control run.4. From these perturbed states, 15 perturbed runs were conducted over several years.5.The evolution of each perturbation was obtained by subtracting the control run from each perturbed run.6.The amplitude of each perturbation was calculated using the energy norm.7. The ensemble average of the amplitudes of these perturbations was obtained.
We compare the evolution of the ensemble average in the UL (70-100 km altitude) and LL (40-70 km altitudes) (Figures 4a and 4b).We observe that the saturation amplitudes of perturbations in the LL (green line) are significantly larger than those in the UL (orange line) (Figure 4a).Furthermore, the former dominates the total amplitude of perturbations, defined within the range of 10-110 km height (blue line).At the onset of the forecasts, since the perturbations are random, their amplitudes quickly decay at the beginning (Figure 4b) because arbitrary random perturbations in the grid domain project mainly on inertia-gravity waves (Toth & Kalnay, 1993).After decaying to some amplitude, they start to grow.The perturbations in the UL experience faster growth (Figure 4b); however, they saturate earlier at a smaller amplitude.This situation is reminiscent of the findings by Toth and Kalnay (1993), who compared the growths of perturbations associated with baroclinic and convective instabilities.In their case, the convective mode is the fast-growing mode, while the baroclinic mode is the slowgrowing mode.However, in our study, the perturbations in the UL and the LL are identified as fast-growing and slow-growing modes, respectively.As explained in the methodology section, larger R amps and longer R ints usually highlight the slow-growing mode (Figure 1a).Therefore, the selection of the R amps and R ints shown in the next subsection is based on the saturation amplitudes of these two modes.
We further compare the average amplitudes of the perturbations in different latitudes in the LL.They exhibit similar trends within a 2-year forecast period (Figure S2 in Supporting Information S1).We obtain the temporally averaged amplitudes of these perturbations from 0005/01/01/06 to 0008/01/01/00 (Figure 4c).In the high latitudes (90 °S-60 °S and 60-90°N), the perturbations exhibit the smallest amplitudes (12%) and are almost symmetric about the equator.In the mid-latitudes (60°S-30°S and 30°N-60°N), the perturbations have larger amplitudes (17%) compared to the high latitudes, and they also exhibit approximate symmetry about the equator.However, the amplitudes of the perturbations in the low latitudes are asymmetric about the equator, with the highest amplitude occurring in 0-30°N (20%).Overall, the amplitudes of the perturbations in the low latitudes are comparable to those in the mid-latitudes.This contrasts with the Earth's atmosphere, where the amplitude of perturbations in the low latitudes (convective mode) is much smaller than that in the mid-latitudes (baroclinic mode) (Toth & Kalnay, 1993).However, we should be cautious about such comparisons.First, since we did not use the nonhydrostatic model nor the convective parameterization scheme, estimating the amplitude of the convective mode is difficult.Second, in this study, we measured the amplitude of the perturbations using the were obtained from the ensemble mean of 15 perturbations from the identical twin experiments, covering the period from 0005/01/01/06 to 0008/01/01/00.The numbers above are the amplitudes (unit: J, the normalization factor is 5*10 17 , and the numbers below are the amplitude percentage at each latitude region to the total amplitude of all regions.Labels A to F correspond to the same labels in Figure 2d. energy norm we defined in Equation 1.However, Toth andKalnay (1993, 1997) used the stream function at 500 hPa as the norm.To our knowledge, we could not find studies using a similar energy norm to ours to measure the saturation and average amplitudes of the convective and baroclinic modes for the Earth's atmosphere.In future studies, we would like to improve the model's representation of the convective mode and use the same norm to compare Venus's and Earth's atmospheres.Finally, notice that the saturation amplitudes of the perturbations in different latitudes in the LL are not significantly different (Figure S2 in Supporting Information S1).As a result, it might be challenging to emphasize some of them by setting different values of R amps and R ints in the breeding cycle experiment.

The Growth Rates and Amplitudes of BVs
Based on the saturation amplitudes of the two modes shown in Figure 4a, different R amps and R ints were used to perform multiple breeding cycles.We analyze the averaged growth rates of BVs from various experiments, considering the period between 0005/02/01/00 and 0005/09/15/18, corresponding to one Venus year (Figure 5).We found that the growth rates of BVs are higher when the R amps are smaller, and the R ints are shorter, aligning with our expectations.Furthermore, we compare the amplitudes of BVs in the LL and the UL while varying the values of R amps and R ints (Figure 5b).When larger R amps (e.g., 100 and 1,000) are used, the amplitudes of BVs in the LL are larger than those in the UL.Conversely, when smaller R amps (e.g., 10) are used, the opposite result is observed.Therefore, we can emphasize different modes by adjusting the values of R amps and R ints .
Next, we examine the time evolutions of the growth rates and amplitudes of BVs in the experiment R amp 1000-R int 06 (R amp = 1,000 and R int = 6 hr).During the initial stages of the breeding cycles (Figure 6a), the growth rates of BV in the UL and LL exhibit significant variations.Because the initial perturbation is random, they decay at the beginning and then start to grow.Subsequently, the growth rate of BV stabilizes, and no significant seasonality is observed, which aligns with our expectations, considering that our model assumes a zero axial tilt for Venus.It would be interesting to investigate the growth rates of BVs after incorporating the actual axial tilt of Venus in future studies, as the planet does possess a slight inclination.In contrast, Mars has a larger axial tilt, resulting in stronger seasonality in the growth rates of BVs within its atmosphere (Greybush et al., 2013).
The amplitude of BV in the LL consistently exceeds that in the UL throughout the entire duration of the R amp 1000-R int 06 experiment (Figure 6b).This phenomenon can be attributed to the R amp value of 1,000, which is significantly larger than the saturation amplitude of the perturbation in the UL and closely comparable to that in the LL (Figure 4a).As the breeding cycle progresses, the perturbations in the UL reach their saturation amplitude and cease further growth, whereas the perturbations in the LL continue to grow.Consequently, the amplitude of the perturbation in the LL eventually dominates the total amplitude of the perturbations.

The Amplitudes and Structures of BVs in Different Regions
To investigate how different instabilities are emphasized by varying the values of R amps and R ints , we selected two experiments for comparison: R amp 10-R int 06 and R amp 1000-R int 06.In this section, our primary focus is on the results obtained from the R amp 1000-R int 06 experiment.The discussion section will cover the findings from the R amp 10-R int 06 experiment.
The latitude-height distributions of the temporally-zonally averaged amplitudes (measured by energy norm) of the BVs reveal that the BV amplitudes in the LL are larger than above and below (Figure 7a).Additionally, we calculated the root mean (temporally-zonally) square of the BV zonal wind (u BV ,m/s), meridional wind (v BV ,m/s), and temperature components (T BV ,K) (Figures 7b-7d).For vertical mass weighting, we multiplied the square root of the difference of sigma levels at each grid point.This normalization approach allows us to visualize the BV components in the lower altitudes, following the conventional methodology used in other studies where the normalization is applied to deviations in control runs (e.g., Sugimoto et al., 2014a;Takagi et al., 2018).However, unlike the calculation of the energy norm of BV, the magnitudes of the BV components are not meridionally weighted by mass, as we did not multiply by the cosine of latitude at each grid point.The magnitudes of the BV components are large in the LL and similar to the deviation fields in the control run (Figure 2).Their distributions are similar as well.Figure 8 illustrates the longitude-height cross-sections of the BV components (u BV ,v BV ,T BV ) at different latitudes.To compare with the deviation in the control, we analyze the BV components at the same time at 0005/06/03/18.In addition, the vertical normalization approach we applied to the deviation in the control is also applied to the BV components, where each grid point value is multiplied by the square root of the difference of sigma levels.In the following paragraphs, we will thoroughly compare the magnitudes and structures of BVs with those of the deviation fields in the control run across different latitudes.Also, Table 1 facilitates the comparison.
In the high latitudes (60-90°N), the amplitude of BV and the magnitudes of BV components at the 40-70 km altitudes are larger than above and below (Figure 7).The magnitudes of BV components (Figures 7b-7d) are comparable to the magnitudes of deviations in the control run (Figures 2a-2c).The longitude-height cross sections of BV (Figures 8a, 8d, and 8g; Movies S1b, S1d, and S1f; Table 1) demonstrate that the phases of u BV and v BV remain nearly constant below 70 km.The phase difference in T BV below and above 65 km is π.u BV and T BV are zonally in phase, while v BV and T BV are zonally out of phase by π/2.These vertical structures are similar to those observed in u′,v′, andT′ in the control run (Figures 3a, 3d, and 3g), although they are not zonally in phase.Furthermore, the latitude-longitude cross-sections of BV at 55 km height (Figure S3 in Supporting Information S1) indicate the presence of a polar vortex structure in high latitude regions (>60°N) in BV.The period of BV is around 7-8 days, similar to the deviation (Table 1).Given the similarity between the structures of BVs and deviations in the control run, we propose that the growth of BVs at the 40-70 km altitudes in the high latitudes is associated with barotropic instability.
In the mid-latitudes (30-60°N), the amplitude of the BVs and the magnitudes of BV zonal and meridional components at 50-70 km are notably larger than above and below (Figures 7a-7c; Movies S2b and S2d).Specifically, Figure 8h shows that the meridional wind and temperature components are slightly out of phase in the zonal direction and titled from lower-east to upper-west at the 50-70 km altitude, similar to those in the deviation
(Table 2).At the 60 km altitude, the temperature deviations and the temperature components in BV are within around ±7 K. Generally, the deviation in the control run and the BV have similar pattern and amplitude (Figure 3), although their phases differ in the zonal direction.Consequently, we can infer that the growth of BV in this region is associated with baroclinic instability.
In the low latitudes (0-30°N), the amplitudes of the BVs at the 52-60 km altitudes exceed those in the mid altitudes (Figure 7a).The patterns in the tropics are connected to those in the mid-latitudes.The longitude-height cross sections reveal that BV at 52-60 km is primarily dominated by wave number one component (Figures 8c,8f,and 8i;Movies S3b,S3d,and S3f), similar to the deviation in the control run (Figures 3c, 3f, and 3i, Table 3), although they are out of phase in the zonal direction.Therefore, we can attribute the growth of BV to the baroclinic instability.The amplitude of BV has a maximum region at the 40-45 km altitude (Figure 7a).The

̅̅̅̅̅ ̅ dσ √
).The BV meridional wind in (d) is superimposed on the BV temperature and shown as contours in (g).(b, e, h) Are similar to (a, d, g) but at 45°N.(c, f, shading in i) Are similar to (a, d, shading in g) but at 0°N.However, the BV zonal wind in (c) is superimposed on the BV temperature (shading in (i)) and shown as contours in (i).magnitude of u BV ,v BV , andT BV at the 40-52 km altitudes are similar, smaller, and larger, respectively, than those at the 52-60 km altitudes (Figures 7b-7d).This pattern resembles the deviations in the control run (Figures 2a-2c).The longitude-height cross sections illustrate that at the 30-52 km altitudes, u BV and T BV are primarily dominated by wave number one component, and they are zonally in phase (Figures 8c,8f,and 8i).Meanwhile, v BV is weak.The phases of u BV and T BV remain constant vertically below 40 km.The structure of BV resembles Kelvin waves in the control run, although they are not zonally in phase (Figures 3c,3f,and 3i).Hence, the growth of BV in this region is attributed to the Rossby-Kelvin instability.
In summary, when R amp is large, the amplitude and structure of BV exhibit similarities to the deviation in the control run, albeit they have different phases in the zonal direction.The growths of BV in the LL in the high latitudes, mid-latitudes, and low latitudes are related to the barotropic, baroclinic, and Rossby-Kelvin instabilities, respectively.Another interesting phenomenon is that the deviation amplitudes in the control run are prominent around the 70 km altitude in the mid-latitudes (Figures 2a-2c).Throughout the evolution of the deviations in the control run (panels a, c, and e in Movies S1-S3), a gradual westward movement of the fields can be observed above the 70 km altitude, which indicates the presence of the thermal tide.However, such distinct maximum regions are not observed in BVs (panels b, d, and f in Movies S1-S3).

Thermal Tide Impact on BVs
Thermal tides commonly occur in the rapidly rotating atmospheres of terrestrial planets such as Earth, Mars, and Venus (Wu et al., 2021).The amplitude of thermal tide in the Earth's atmosphere is small (Guerlet et al., 2023).For the thermal tide in the Martian atmosphere, dust and water ice are important excitation sources because they absorb the visible sunlight and the infrared radiation emitted from the ground (Wu et al., 2020).In the Venus atmosphere, solar heating in the thick cloud layer excites the thermal tide (Apt et al., 1980;Kouyama et al., 2019).Since thermal tide is important for the general circulation of the Venus atmosphere (Takagi et al., 2018), it would be interesting to investigate how BVs develop under its influence.We analyze the composite means of deviations in the control run and the magnitudes of BV components.The composite means are obtained following the subsolar position .Similar plots for the Qz case (thermal tide is excluded) are presented in Figures S4-S7 in Supporting Information S1 to provide a basis for comparison.The central position (180°) in the composite mean plot represents the subsolar position.Given that the Venus rotation is set counterclockwise in this study, 0-180°corresponds to the morning hemisphere, while 180-360°corresponds to the evening hemisphere.The composite means of the deviations in the control run (panels a, c, and e in Figures 9-12) exhibit similarities to the findings reported by Takagi et al. (2018).
The longitude-height cross sections at the 45°N latitude of the composite means show that the magnitudes of the BV components are greater in the morning hemisphere at the 60-70 km altitudes (Figures 9b, 9d, and 9f).Correspondingly, as shown in the latitude-longitude distributions of BV magnitudes at the 65 km altitudes of the composite means, the local maximum magnitudes in the morning hemisphere in 30-60°N extend from northeast to southwest (Figure 12).At the 60 km altitude (Figure 11), the local maximum magnitudes of BV zonal wind and temperature components can also be observed on the morning side.In contrast, the equivalent plots for the Qz case demonstrate zonally uniform magnitudes (Figures S4-S7 in Supporting Information S1).Therefore, the nonuniform zonal distribution in the Qt case suggests that the thermal tide increases the growth of BV in the morning hemisphere in the mid-latitudes.
To explain this phenomenon, we analyze the composite means of the temperature deviations in the control run (Figures 11e and 12e).They indicate that thermal tide results in higher temperatures in the low latitudes and lower temperatures in the mid-latitudes in the morning hemisphere around 0-90°E (Figures 11e and 12e).The temperature distribution influenced by the thermal tide may increase the mean available potential energy in the midlatitudes and subsequently enhance the baroclinic instability.Consequently, the amplitude of BV is greater in this region and extends downstream to the east (90-180°E).
In comparison, the longitude-height cross sections at around 0°latitude of the composite means show that the magnitudes of the BV components are relatively uniform in the zonal direction at the 40-60 km altitudes (Figures 10b, 10d, and 10f), suggesting a minimal influence of the thermal tide on BV growth in this region.Above the 60 km altitude, the thermal tide has wave number two (semi-diurnal tide, Takagi et al., 2018)) at low latitudes (Figures 10a, 10c, and 10e).However, in the same region, as we analyzed before, the breeding cycle in this experiment (R amp 1000-R int 06) could not generate a significant BV (Figure 7 and Table 3).Therefore, we Longitude-height distributions of the squared values of (b) BV zonal wind (dσu BV 2 comp , m 2 s 2 ), (d) BV meridional wind (dσv BV 2 comp , m 2 s 2 ), and (f) BV temperature (dσT BV 2 comp , K 2 ) components associated with the thermal tides at 45°N from the R amp 1000-R int 06 experiment, weighted by the difference of the sigma level dσ.All distributions are calculated using composite means over 1 Venus year (time 0005/02/01/00 to 0005/09/15/18) at the subsolar point (fixed at the center of each panel).Since the planet rotates from left to right in the present study, 0-180°and 180-360°represent the morning and evening hemispheres, respectively.(u BV 2 comp , m 2 s 2 ), (d) BV meridional wind (v BV 2 comp , m 2 s 2 ), and (f) BV temperature (T BV 2 comp , K 2 ) components associated with the thermal tides at 60 km altitude from the R amp 1000-R int 06 experiment.All distributions are calculated using composite means over 1 Venus year (time 0005/02/01/00 to 0005/09/15/18) at the subsolar point (fixed at the center of each panel).Since the planet rotates from left to right in the present study, 0-180°and 180-360°represent the morning and evening hemispheres, respectively.We propose that these small-scale features might be related to the generation and breaking of the gravity waves.Sugimoto et al. (2021) suggested that in the Venus atmosphere, gravity waves above 70 km in the low, mid, and high latitudes could be generated from the thermal tide, baroclinic waves, and barotropic waves, respectively, through the spontaneous mechanism (Yasuda et al., 2015).It is plausible that the small-scale features above the 70 km altitude in the BV might be the gravity waves generated by a similar mechanism.As they are generated, energy is transferred from vortices to gravity waves, leading to the growth of these smaller-scale features.Furthermore, as the gravity waves propagate to higher altitudes with lower density, their amplitudes tend to increase.These larger amplitude gravity waves could create an unstable layer, promoting conditions of convective instability.This convective instability could then facilitate the rapid growth of BV, indicating that the gravity waves eventually break, although convective instability is treated in a very simple way by the convection adjustment scheme in the model.
Composite mean plots of the squared values of the BV zonal wind and temperature components indicate that there are two bands of maximum values in the mid-latitudes and tropics at the 80 km altitude in the R amp 10-R int 06 experiment (Figures S8b and S8f in Supporting Information S1).These regions closely align with the jet exit areas observed in the composite mean plots of the zonal wind deviation in the control run (from the left to the right of Figure S8a in Supporting Information S1, the regions where the zonal wind deviations change from positive to negative).The longitude-height cross sections at 0°and 45°N are shown in Figure S9 in Supporting Information S1.Two maximum magnitudes of BV zonal wind and temperature components can also be seen in the midlatitudes and tropics at an altitude above 70 km (Figures S9a,S9b,S9e,and S9f in Supporting Information S1).Based on previous studies (e.g., Sugimoto et al., 2021), it is known that gravity waves can be generated by the thermal tide at the jet exit regions through a spontaneous mechanism.Notably, such bands of maximum values are absent in the Qz case (Figure not shown).Therefore, composite means of the BVs in the Qt case suggest a potential relationship between the gravity waves and the thermal tide.
Baroclinic waves can generate gravity waves above 70 km altitude in the mid-latitudes through a spontaneous mechanism.Given the stronger baroclinic instability below 70 km in the morning hemisphere in the mid-latitudes observed in the composite mean plots (Figures 9b,9d,and 9f), it is reasonable to observe that the generation of gravity waves around 70 km altitude and their breaking above 80 km altitude are more pronounced in the morning hemisphere in the mid-latitudes (Figures S8b,S8d,S8f,S9a,S9c,and S9e in Supporting Information S1).These findings provide valuable insights into the generation and breaking of gravity waves in both mid-latitudes and tropical regions.

Predictability
In this section, we describe the error growth, doubling (e-folding) time, and predictability in the mid-latitudes (30-60°N) at around the 54-km altitude.Since our model is hydrostatic, it did not explicitly resolve the convection.According to Zhang et al. (2003), moist convection is the key process that drives rapid initial error growth for the Earth atmosphere.In general, the error doubling time will be shorter if the model is more realistic (Judt, 2018).Therefore, we expect that the error growth, doubling (e-folding) time, and predictability in this study could change if a more realistic model, such as a nonhydrostatic one, is developed for the Venus atmosphere.Despite the limitation, we conducted additional experiments to study this topic after finishing the breeding cycles.
From the R amp 10-R int 06 experiment, we randomly picked perturbed forecasts at 59 different times and extended them to several months.
The kinetic energy (KE) and the total energy (TE) were defined in Equations 3 and 4, following similar metrics in Judt (2018).
The meaning of the variables is the same as those described in the methodology.The unit of KE is m 2 s 2 and the unit of TE is J/kg.For each forecast, their values at each grid point at the 54-km altitude in 30-60°N were averaged, weighted by the mass of each grid point.These time series were then averaged to obtain the average time series for KE and TE (Figure 14 for KE.TE is similar but not shown).
Figure 14 shows that the error growth rates differ at different lead times.This time series can be approximately divided into three sections: days 0-17 (slow), days 17-32 (fast), and days 32-60 (slow).On days 0-17, the perturbation has small-scale features (not shown).We suggest that they are related to the radiation of some gravity waves.Since the perturbation we gave will contain mostly unbalanced gravity waves, the balance component that grows by baroclinic instability will be small, so it takes around 17 days to be larger than the gravity waves.On days 17-32, the growth rate is higher, and the corresponding average doubling time and e-folding time are around 4.6 and 6.6 days, respectively.During this period, the perturbation amplitude could be mainly developed by the baroclinic instability.In comparison, the efolding time of the kinetic energy for wave numbers 5 and 6 at 30-60°N at 54 km altitude is around 5.4 days in Sugimoto et al. (2014a) (their Figure 8), comparable to our number.After day 32, the growth rate is slower, indicating it is approaching the saturation level.
Predictability includes "intrinsic" predictability and "practical" predictability (Lorenz, 1996;Sun & Zhang, 2016).The "intrinsic" refers to the limit of prediction if the initial condition and the forecast models are almost perfect.In contrast, the "practical" refers to the limit of prediction using the initial condition derived from the current optimal analysis along with the best available forecast model.Given that the initial perturbation (∼0.1 m 2 s 2 ) we assigned is not "very small" (<10 5 m 2 s 2 ), we infer that the "intrinsic" predictability might be longer than 1 month (Figure 14).However, the typical analysis and observational errors would be much larger than this initial perturbation, given the very limited number of observations for the Venus atmosphere.In addition, our model is not perfect.Therefore, the "practical" predictability may be much shorter than the "intrinsic" predictability.
We compare the situations between the Venus and Earth's atmospheres.The error doubling time for the Earth's atmosphere is around 1.2-1.7 days in the last two decades (Judt, 2018).Using a global convection-permitting model, Judt (2018) found that the average doubling times of the total energy over the core of the baroclinic phase in the Earth's atmosphere (0-11 km) is around 1.6 days (their Figures 8 and 11).Compared to the 4.6 days error doubling time in our Venus GCM, it indicates that the development of the baroclinic instability in the Venus atmosphere may be slower than that of the Earth's atmosphere.The "intrinsic" predictability of the Earth's tropospheric atmosphere is around 2-3 weeks (Judt, 2018), which might be shorter than our Venus GCM.Again, we should be cautious about such comparisons because the region used to calculate the energy is different, and more realistic models might have different results.

Challenges
In Section 3, we first discussed the BVs from the experiment R amp 1000-R int 06, and compared them to the deviations from the zonal mean in the control run.Since R amp = 1,000 emphasizes the slow-growing modes in the LL (40-70 km altitude).We found that the structure and amplitude of BVs resemble those of the deviations.Such results are reasonable because the waves in the control run and the BVs grow due to the same instabilities: barotropic, baroclinic, and Rossby-Kelvin instabilities.We then analyzed the BVs from the experiment R amp 10-R int 06.Because R amp = 10 emphasizes the fast-growing modes in the UL (70-100 km altitude), we found smallscale features that could be related to the generation of gravity waves through the spontaneous mechanism and the breaking of gravity waves due to convective instability.In short, BV can identify the growing mode due to instability, spontaneous emission, and breaking of gravity waves.
The thermal tide can be seen in the control run by taking the composite mean of the deviations following the sun position, but it cannot be seen in BV because the thermal tide is a neutral wave that does not grow.However, the thermal tide has an indirect impact on BV.Firstly, it may intensify the growth of BV on the day side at midlatitudes in the LL.Secondly, thermal tide generates gravity waves, which BV can detect.One of the limitations of our study is that, since our model is hydrostatic, it does not resolve the convection explicitly.Although the convective adjustment is implemented, it may be difficult to represent the convective instability realistically.If using the nonhydrostatic model, we expect that BV can identify the convective instability in the cloud layer as well, and we can compare them to those in the Earth's atmosphere.

Conclusion
Previous studies using the AFES-Venus model have investigated the barotropic, baroclinic, and Rossby-Kelvin instabilities in the Venus atmosphere.In this study, we utilized the breeding of growing modes (BGM) method to enhance our understanding of the instabilities in the Venus atmosphere.We first conducted the control run using the AFES-Venus model to simulate the Venus atmosphere.Although we did not explicitly separate the fast and slow waves in the analysis, we found that the basic state and the deviations from the zonal mean in the control run are similar to the previous studies that used similar model settings and mainly focused on fast waves.We also observed the barotropic, baroclinic, and Rossby-Kelvin instabilities in our simulations.
We conducted identical twin experiments to investigate the saturation amplitudes of the fast-growing and slowgrowing modes.We added 15 random perturbations to the initial condition of the control run to initiate 15 perturbed run experiments.The ensemble average of the amplitudes of these perturbations revealed that, in the beginning, the perturbations in the upper layer (UL, 70-100 km altitudes) grow faster compared to those in the lower layer (LL, 40-70 km altitudes).However, the former saturates early at a significantly low amplitude, indicating a fast-growing mode, while the latter represents a slow-growing mode.By analyzing the temporal average of the perturbations in the LL across different latitudes, we found that the amplitudes of perturbations in the high latitudes are the smallest.The amplitudes in the mid and high latitudes are symmetric about the equator, while those at the low latitudes are asymmetric.The amplitudes in the low latitudes are comparable to those in the mid-latitudes.This result contrasts with the Earth's atmosphere (Toth & Kalnay, 1993), where the amplitude of the convective mode in the low latitudes is much smaller than the amplitude of the baroclinic mode in the mid-latitudes.However, we should be cautious about such comparisons because our model may not represent the convection very well, and we measured the amplitude of the perturbations using a different norm.
Based on the saturated amplitudes of the fast and slow-growing modes observed in the identical twin experiments, we set different rescaling amplitudes (R amps ) and rescaling intervals (R ints ) to conduct multiple breeding cycle experiments.We found that the growth rates of BVs show small seasonal variations.The growth rates are higher when smaller R amps and R ints are used.When larger R amps are used, the amplitudes of the BVs in the LL are comparable to the amplitudes of the deviations in the control run, and they are larger than those in the UL.Moreover, the structures of the BV zonal wind, meridional wind, and temperature components in the LL resemble the deviations in the control run despite having zonally different phases.This similarity can be attributed to the growths of deviations in the control run and BV due to the barotropic, baroclinic, and Rossby-Kelvin instabilities in the high, mid, and low latitudes, respectively.
The thermal tide has an indirect impact on the growth of BVs.By analyzing the composite means following the sub-solar position, it was observed that at the 60-70 km altitudes in the mid-latitudes, the magnitudes of the BV components are larger in the morning hemisphere compared to those in the evening hemisphere.This finding suggests that the thermal tide intensifies the baroclinic instability and the growth of BVs in the morning hemisphere in this region.
When using smaller R amps , we observed the presence of small-scale features above the 70 km altitude.Based on our analysis, we propose that these small-scale features may indicate the generation and breaking of gravity waves.
Based on the extended forecasts from the perturbed forecasts during the breeding cycle, the error doubling time and predictability related to the baroclinic instability in the mid-latitudes are investigated.Measured by the kinetic energy, the error doubling time is around 4.6 days, much longer than in Earth's atmosphere.The intrinsic predictability would be longer than 1 month.
As highlighted in the introduction, BV helps elucidate the physical origins of the instabilities.BV energy equations have been used to study instabilities in various systems, such as the tropical Pacific Ocean and the Martian atmosphere.Furthermore, by utilizing multiple BVs obtained from breeding cycles, the BV dimension (BVD) can be calculated, providing a measure of the effective local dimensionality of a chaotic system.Drawing from the insights gained in this study, we will extend our understanding of the instabilities in the Venus atmosphere using BV energy equations and BVD analysis in future studies.

Figure 1 .
Figure 1.(a) Amplitude evolutions of the slow-growing mode (black solid line) and the fast-growing mode (red solid line).The dashed lines represent the rescaling amplitudes (R amps ).This figure is adapted from Figure 7 in Toth and Kalnay (1993).(b) Steps involved in conducting breeding cycle experiments to obtain bred vectors (BVs).

Figure 2 .
Figure 2. Latitude-height distributions of the variables in the control run.(a) Zonally and temporally averaged zonal wind (contour, ms 1 ) and the RMS (square root of the zonally-temporally averaged square) of zonal wind deviations from the zonally averaged zonal wind ( ̅̅̅̅̅ ̅ dσ √ ̅̅̅̅̅̅ ̅ u′ 2 √ ,ms 1 ), weighted by the square root of the difference of the sigma level ( ̅̅̅̅̅ ̅ dσ √ ) (color shade, ms 1 ).(b) Is similar to (a) but for meridional wind (ms 1 ).(c) The color shades are similar to those in (a) but for temperature (K).The zonally and temporally averaged temperature was first obtained, and the deviations from the meridional average were shown in contours (solid lines indicate positive values while dashed lines indicate negative values).The temporal average covers one Venus year from 0005/02/01/00 to 0005/09/15/18.(d) Indicates different regions.UL represents the upper layer at 70-100 km altitudes, and LL represents the lower layer at 40-70 km.In LL, it is divided into 6 regions with a meridional interval of 30°.Labels A to F correspond to the same labels in Figure 4c.

Figure 3 .
Figure 3. Longitude-height distributions of the (a) zonal wind ( ̅̅̅̅̅ ̅ dσ √ u′, ms 1 ), (d) meridional wind ( ̅̅̅̅̅ ̅ dσ √ v′, ms 1 ), and (g) temperature ( ̅̅̅̅̅ ̅ dσ √ T′, K, shading) deviations , shading), the amplitudes of u′ and v′ are relatively large and have connections with the mid-latitudes.The longitude-height cross sections in the low latitudes(Figures 3c, 3f, and 3i; Movies S3a, S3c, and S3e)  show that the deviation fields in the control run at the 52-60 km altitude are primarily dominated by wave number one component.AsTakagi et al. (2022) suggested, the 5.8-day waves above the critical level (∼55 km in their study) in the low latitudes are associated with the midlatitude Rossby waves.Similarly, the waves identified in our study at the 52-60 km altitude are also associated with the mid-latitude Rossby waves related to the baroclinic instability.At the altitudes of 40-52 km, the amplitudes of u′, v′, and T′ are similar, smaller, and larger respectively, in comparison to those at the 52-60 km altitudes (Figures2a-2cshading).The longitude-height cross sections(Figures 3c, 3f, and 3i; Movies S3a, S3c,  and S3e) demonstrate that at the 30-52 km altitude, u′ and T′ are dominated by wave number one component and they are zonally in phase.The altitude of T′ is slightly lower than u′ for the same zonal phase.v′ is weaker than those above.Furthermore, the phases of u′ and T′ remain constant vertically below the 40 km altitude, indicating the barotropic nature of the waves.In comparison, such altitude is around 42 km inTakagi et al. (2022).The wave characteristics resemble the Kelvin waves observed inTakagi et al. (2022) (Table

Figure 4 .
Figure 4. (a) Time evolution of perturbation amplitudes (unit: J, the normalization factor is 5*10 17 ) at different vertical ranges: UL (70-100 km) (orange), LL (40-70 km) (green), and 10-110 km (blue), calculated from the ensemble mean of 15 perturbations from identical twin experiments.The x-axis indicates forecast days since 0005/01/01/00.(b) Is similar to (a), but only the first 90 days are shown, and the y-axis is in the log scale.(c) Average perturbation amplitudes at different latitudes in LLwere obtained from the ensemble mean of 15 perturbations from the identical twin experiments, covering the period from 0005/01/01/06 to 0008/01/01/00.The numbers above are the amplitudes (unit: J, the normalization factor is 5*10 17 , and the numbers below are the amplitude percentage at each latitude region to the total amplitude of all regions.Labels A to F correspond to the same labels in Figure2d.

Figure 5 .
Figure 5. (a) The temporal average of the daily growth rates of the BVs obtained from various experiments with different R amps (unit: J, the normalization factor is 5*10 17 ) and R ints (hour).(b) The temporal average of the amplitudes (unit: J, the normalization factor is 5*10 17 ) of the BVs in LL (40-70 km) is subtracted from those in the UL (70-100 km), obtained from different experiments with different R amps and R ints .The average period is 1 Venus year from 0005/02/01/00 to 0005/09/ 15/18.

Figure 6 .
Figure 6.(a) Time series of the daily growth rates of BV in LL (40-70 km) (orange) and UL (70-100 km) (blue) from the experiment R amp 1000-R int 06.The dots represent the 6-hourly data, while the lines depict the ±5-day moving averages.The xaxis indicates days since 0005/01/01/00.(b) Similar representation as (a), but for the amplitudes (unit: J, the normalization factor 5*10 17 ) of the BVs shown every 6 hr.

Figure 11 .
Figure 11.Longitude-latitude distributions of (a) zonal wind deviations (u′ comp , ms 1 ), (c) meridional wind deviations (v′ comp , ms 1 ), and (e) temperature deviations (T′ comp , K) associated with the thermal tides at 60 km altitude from the control run.Longitude-latitude distributions of the squared values of (b) BV zonal wind

Figure 14 .
Figure14.Average amplitudes of 59 perturbations (solid black line) (unit: J, the normalization factor is 5*10 17 ).The amplitude is measured by the massweighted average of the kinetic energy at 54 km altitude in 30-60°N.The 59 forecasts are the extension of the perturbed forecasts at random times during the breeding cycle from the experiment R amp 10-R int 06.The X-axis indicates the forecast days.The gray color indicates 1 standard deviation (from 59 perturbations) above and below the average.The y-axis is in the log scale.

Table 1
Comparison Between the Deviations and BVs for the High Latitude Regions

Barotropic instability, cross-section at 80°N
Note.The BVs are from R amp 1000-R int 06 experiment.U, V, and T indicate zonal wind, meridional wind, and temperature, respectively.

Table 2
Similar to Table 1 but for the Mid-Latitude Regions

Table 3
Similar to Table 1 but for the Low-Latitude Regions