Single-particle Lagrangian statistics from direct numerical simulations of rotating-stratified turbulence

Geophysical fluid flows are predominantly turbulent and often strongly affected by the Earth's rotation, as well as by stable density stratification. Using direct numerical simulations of forced Boussinesq equations, we study the influence of these effects on the motion of fluid particles, focusing on cases where the frequencies associated with rotation and stratification (RaS), N and f respectively, are held at a fixed ratio $N/f=5$. As the intensity of RaS increases, a sharp transition is observed between a regime dominated by eddies to a regime dominated by waves, which can also be seemingly described by simply comparing the time scale 1/N and $\tau_\eta$ (the Kolmogorov time scale). We perform a detailed study of Lagrangian statistics of acceleration, velocity and related quantities in the two regimes. The flow anisotropy induces a clear difference between particle motion in the horizontal and vertical directions. In the regime $N\tau_\eta<1$, acceleration statistics in both horizontal and vertical directions, exhibit well known characteristics of isotropic turbulence. In contrast for $N\tau_\eta>1$, they are directly influenced by imposed RaS. The Lagrangian velocity statistics exhibit visible anisotropy for all runs; nevertheless the degree of anisotropy becomes very strong in the regime $N\tau_\eta>1$. We find that in the regime $N\tau_\eta<1$, rotation enhances the mean displacement of particles in horizontal planes at short times, but inhibits them at longer times. This inhibition of horizontal displacement becomes stronger for $N\tau_\eta>1$, with no clear diffusive behavior. Displacements in the vertical direction are always inhibited. The inhibition becomes extremely strong when $N\tau_\eta>1$, with the particles almost being trapped horizontally.

(Dated: September 30, 2019) Geophysical fluid flows are predominantly turbulent and often strongly affected by the Earth's rotation, as well as by stable density stratification. Using direct numerical simulations of forced Boussinesq equations, we study the influence of these effects on the motion of fluid particles, focusing on cases where the frequencies associated with stratification and rotation, N and f respectively, are held at a fixed ratio N/f = 5. The simulations are performed in a periodic domain, at Reynolds number Re ≈ 4000, and Froude number F r in the range 0.03 F r 0.2 (with Rossby number Ro = 5F r). As the intensity of rotation and stratification increases, a sharp transition is observed between a regime dominated by eddies to a regime dominated by waves, which corresponds to F r 0.07. For the given runs, this transition to a wave-dominated regime can also be seemingly described by simply comparing the time scales 1/N and τ η , the latter being the Kolmogorov time scale based on the kinetic energy dissipation. We perform a detailed study of Lagrangian statistics of acceleration, velocity and related quantities in the two regimes. The flow anisotropy induces a clear difference between particle motion in the horizontal and vertical directions. In the regime N τ η < 1, acceleration statistics in both horizontal and vertical directions, exhibit well known characteristics of isotropic turbulence, such as probability density functions (PDFs) with wide tails and acceleration variance approximately scaling with Kolmogorov scales. In contrast for N τ η > 1, they behave very differently, exhibiting direct influence of the imposed rotation and stratification. On the other hand, the Lagrangian velocity statistics exhibit visible anisotropy for all runs; nevertheless the degree of anisotropy becomes very strong in the regime N τ η > 1. In considering the mean displacement of particles, we find that in the regime N τ η < 1, rotation enhances the displacements in horizontal planes at short times, but inhibits them at longer times, however a transition to diffusive behavior is evident for all cases. This inhibition of the horizontal displacements becomes stronger in the regime N τ η > 1, with no clear diffusive behavior. On the other hand, displacements in the vertical direction are always inhibited, but the inhibition becomes extremely strong in the N τ η > 1 regime, leading to a scenario where particles almost appear to be trapped in horizontal planes.

I. INTRODUCTION
Transport of material substances plays a crucial role in many geophysical processes [1,2], e.g. dispersion of pollutants and contaminants [3,4], droplet dynamics in clouds [5], mixing of planktons and other biomatter in the oceans [6]. The Lagrangian viewpoint following the motion of fluid particles [7] or analogous entities such as Brownian particles or inertial particles [8,9], provides a natural description of such transport processes. Not only are most geophysical flows turbulent, they are also often strongly influenced by anisotropies due to effects such as rotation (Coriolis force) and stratification (buoyancy force). While particle dispersion has been studied extensively in isotropic turbulence, e.g. see reviews [10][11][12], it has only recently started receiving attention in anisotropic flows. In particular, several studies have focused on flows considering either the effects of rotation [13][14][15][16] or stratification separately [17][18][19]. However, in many applications, adequately describing the observed flow physics necessitates examining the combined effects of rotation and stratification (RaS), e.g. in the southern abyssal oceans with particularly high mixing intensities [20][21][22].
Arguably the most challenging aspect of studying the combined effects of RaS in turbulence is the enormous range of spatial and temporal scales associated with such flows. In isotropic turbulence, typically the only governing parameter is the Reynolds number (Re), which directly provides a measure of the range of scales in the flow. However, at least three additional parameters have to be considered in rotating-stratified turbulence, namely, the Rossby number (Ro) and the Froude number (F r) which respectively measure the strength of RaS, and the Prandtl number (P r) which is the ratio of the fluid viscosity to the thermal diffusivity (see Section II for precise definitions). This leads to consideration of additional length and time scales -beyond the already wide range existing in isotropic turbulence -and makes the flow dynamics far more involved. In particular, it has been recognized that the interaction between linear and non-linear processes in such flows leads to a rich variety of complex behavior, such as spontaneous generation of helicity, dual cascading of energy and many more [23][24][25]. Advances in computing power in the last decade, has allowed for significant progress in understanding such flows, although mostly from the Eulerian perspective [26][27][28][29][30][31].
In this work, utilizing direct numerical simulations (DNS) of Boussinesq equations, our objective is to investigate the dispersion of fluid particles in rotating-stratified turbulence. In particular, we focus on the motion of single individual particles. In view of the very large parameter space associated with rotating-stratified turbulence, we limit our investigation to the case where Ro/F r = 5, i.e., the strength of stratification is five times that of rotation, which is largely relevant to the southern abyssal oceans [20,21]. We use the same grid sizes for all cases, giving Re ≈ 4000 − 5000, and Froude numbers in the range 0.03 F r 0.2. This regime of F r corresponds to the transitional regime, where both turbulent eddies and inertia-gravity waves are expected to play an important role [25]. As recently demonstrated in the case of purely stratified flows, the transition towards a wave dominated regime is accompanied by intermittent large-scale vertical drafts [32,33]. We show that the presence of additional rotation does not, in general, prevent the manifestation of this intermittency and thereby also investigate its effect on Lagrangian statistics.
In studying the motion of single particles, we investigate the properties of both their acceleration and velocity. The intrinsic anisotropy of the flow makes it necessary to distinguish throughout between motion in the horizontal plane and in the vertical direction, which are respectively perpendicular and parallel to the direction of imposed stratification (and also the axis of rotation). For acceleration statistics, which are reflective of small-scales of turbulence, two distinct regimes are observed depending on the strength of imposed RaS. At weak or moderate RaS, the properties of acceleration are qualitatively similar to those documented in isotropic turbulence (with minor quantitative deviations), suggesting that the imposed anisotropy at large scales does not signifi-cantly affect the small-scales. However, when RaS becomes very strong, acceleration statistics are strongly modified. This is reflected, among others, by a very strong anisotropy between the horizontal and vertical motions. In the range of parameters considered here, we observe a sharp transition between the two regimes, which can be simply quantified with the ratio of the Kolmogorov time scale and the stratification time scale (which is simply the inverse of the Brunt-Väisäla frequency). The qualitative change in the statistical properties of acceleration occurs simultaneously with the appearance of intermittent bursts in the flow.
The Lagrangian velocity statistics, in contrast to acceleration, are always affected by the imposed RaS, with the degree of anisotropy expectedly increasing with strength of imposed RaS. However, the transition observed for acceleration statistics is also visible in velocity statistics, resulting in even stronger anisotropy in velocity statistics when RaS is very strong. Accordingly, we characterize the integral time scales based on velocity autocorrelations and investigate the displacement of particles. The particles move ballistically at short times. At long times, the emergence of a diffusive regime is evident for runs with moderate RaS. However for strong RaS, the dispersion, particularly in vertical direction is strongly suppressed, qualitatively consistent with earlier works on purely stratified flows.
This paper is organized as follows. In Section II, we discuss our numerical methods. Essential features of the underlying Eulerian flow are discussed in Section III. The statistics of acceleration are presented in Section IV. Section V presents our results on the Lagrangian autocorrelation functions of velocity along with results on how particles are displaced by the flow. Finally, Section VI contains our concluding remarks.

II. NUMERICAL METHOD AND DATABASE
We simulate the Eulerian flow by numerically integrating the incompressible Boussinesq equations in a rotating frame, with constant solid body rotation rate Ω (and frequency f = 2Ω) and gravity g anti-aligned in the vertical (z) direction: Here, u = (u, v, w) is the velocity field, θ is the temperature fluctuation (in units of velocity), P is the pressure normalized by background density, ν is the kinematic viscosity, κ is the thermal diffusivity and N is the Brunt-Väisäla frequency, which characterizes the strength of imposed stratification. We take ν = κ, assuming the Prandtl number to be unity. A stochastic forcing term, Φ, is utilized to achieve and maintain a statistically stationary state. The forcing is random in time and isotropic in Fourier space, with the energy being injected at large scales in a spherical shell of wavenumbers given by 2 ≤ |k| ≤ 3 [34].
To characterize rotation and stratification, we introduce the dimensionless parameters, the Rossby (Ro) and Froude (F r) numbers defined as: where U and L are respectively the characteristic velocity and length scale of the large-eddies of the flow. An important parameter is the ratio N/f (= Ro/F r), which measures the relative strength of rotation and stratification. As already mentioned, we maintain N/f = 5 for all runs, as relevant in some oceanographic situations [20,21].
In addition, we define the following dimensionless numbers where ǫ is the mean dissipation rate of turbulent kinetic energy and ω 2 is the mean enstrophy density (with ω being the magnitude of vorticity). It is worth noting that in homogeneous turbulence, as will be considered in this work, the two are related as ω 2 = ǫ /ν = 1/τ 2 η , where is the Kolmogorov time scale characterizing the motion of turbulent eddies [35]. The Reynolds number Re is defined in the classical way. The parameter R IB , often referred to as the buoyancy Reynolds number, provides a direct measure of the relative importance of the waves induced by stratification with respect to turbulence in the flow. This is evident by considering the Ozmidov length scale ℓ OZ = ( ǫ /N 3 ) 1/2 , which characterizes the gravity waves and the Kolmogorov length scale, η = (ν 3 / ǫ ) 1/4 , which characterizes the dissipation scale of turbulence, thereby giving R IB = (ℓ OZ /η) 4/3 . Alternatively, R IB can be written as the ratio of two time scales: where T N = 1/N is the stratification time scale (and τ η is the Kolmogorov time scale), once again characterizing the relative strength of stratification to that of turbulent eddies. An alternative definition of the buoyancy Reynolds number [25,28] is given by the parameter R B . The two parameters are related as R IB = βR B , where β = ǫ L/U 3 . The ratio β is expected to be constant as a consequence of the dissipation anomaly in turbulent flows [36]. In practice, β appears to depend on the strength of stratification, especially when the Reynolds number is not very high [25]. In our analysis, we will utilize R IB to quantify the relative strength of stratification and turbulence. The relevance of this choice will become evident as our results are presented in later sections.
Finally, in addition to the Rossby number Ro based on large scales, it is useful to define the micro-Rossby number Ro ω , which in contrast to R IB , measures the relative strength of small-scale turbulent motions to that of the imposed rotation. However, using ǫ = ν ω 2 due to statistical homogeneity, Ro ω can be simply related to R IB as: Since N/f is held constant in this work, the parameters Ro ω and R IB essentially provide the same information.
The database utilized in the current work, along with the main simulation parameters, is summarized in Table I. The simulations were carried out using the Geophysical High-Order Suite for Turbulence (GHOST) code, a versatile, highly parallelized, pseudo-spectral code, utilizing hybrid MPI-OpenMP programming model [37], with second order explicit Runge-Kutta time stepping. All runs correspond to a (2π) 3 periodic domain with 512 3 grid points. The parameters N and f are varied over a range of values keeping the ratio N/f fixed at 5. In addition, for a systematic comparison, we have also performed additional simulations by setting N = 0 or f = 0 or both. The simulation with N = f = 0 (case 0) simply corresponds to homogeneous isotropic turbulence (HIT). In estimating the various dimensionless numbers defined earlier, we have used L = L f , where L f = 2π/(2.5) corresponds to the kinetic energy injection scale, and U = |u| 2 1/2 , the mean amplitude of the velocity fluctuations.
The runs discussed here were started from an initial condition, consisting of a few random modes in the velocity field, whereas the temperature field, θ, was initialized to zero. The Boussinesq equations were integrated until a statistically stationary state was reached. Fig. 1 shows the values of the energy dissipation, averaged over space, and as a function of time. As typical of numerical forced simulations, we observe that the mean dissipation rate starts from zero and rapidly increases to a large peak value, until finally reaching a (plausibly) statistically stationary state around a mean value lower than the peak. The particles are injected at the onset of the stationary state, as indicated by a cross symbol in Fig. 1 period of at least 10 − 20T E in the stationary state for each case, where T E = L/U is the large eddy turnover time.  Table I. Time t has been rescaled by the large eddy turnover time T E = L/U . Fluid particles were injected at the time indicated by the crosses. In Table I, an important point to note is that the values of R B and R IB monotonically decrease as N and f are increased. Values of N and f higher than those shown are avoided, since the resulting regime would be dominated by waves [25]. The simulations with either N = 0 or f = 0 have their f or N value (respectively) corresponding to rotating-stratified run 5, with one of the largest N and f values. However, as evident from the table the values of the F r (and Ro) numbers corresponding to runs 5 and 7 (respectively runs 5 and 8) differ significantly, despite identical values of the stratification (respectively rotation rate). These variations are a consequence of our dynamical definition of U , which accounts for the subtle interplay between turbulence, rotation and stratification. In particular, they reflect the significant variation in U itself, which will be discussed in Section III. Importantly, we observe that the value of R IB is lower than 0.5 for Runs 4-6, as well as for the purely stratified run, 7, suggesting a dominant role of the waves in these runs. In contrast, R IB > 2 in all other runs, possibly pointing to a prevalent role of the eddies. This leads to two qualitatively different behaviors, as demonstrated by our results in the following sections.
In order to obtain the Lagrangian statistics, the fluid particles are tracked in time along with the Eulerian flow, according to the basic equation of motion where the superscript + denotes a Lagrangian quantity and the fluid particle velocity is simply defined to be the Eulerian velocity field evaluated at the instantaneous particle position. For each run listed in Table I, we additionally tracked the motion of 1.5M particles (M=10 6 ), which are randomly distributed throughout the entire domain at the time of injection. The number of fluid particles are held constant for all cases, since given the same grid size and approximately similar Reynolds numbers, their sampling requirements are also approximately the same. Similar to the Eulerian grid, the particles are distributed among parallel processors and tracked in time together with the velocity field using a second order Runge-Kutta scheme, and interpolation methods based on cubic splines [38][39][40]. As mentioned earlier and marked in Fig. 1, the particles were injected after the Eulerian flow reached a statistically stationary state.

III. EULERIAN VELOCITY FIELD
In this section, we provide a brief overview of the underlying Eulerian flow, which establishes the necessary framework for better understanding the Lagrangian statistics reported subsequently. We begin by considering the energy spectra, as shown in Fig. 2. The spectra, normalized using the large scale variables U and L, are split into two sub-figures depending on the strength of stratification. Fig. 2a shows the spectra for the runs with relatively weak rotation and stratification (RaS), corresponding to F r ≥ 0.086 or equivalently R IB > 1 (including the HIT and rotation-only runs). The main observation is that all the spectra demonstrate approximate k −5/3 scaling in an intermediate range of k, following Kolmogorov's theory [41]. This is expectedly most evident for the HIT run and the rot-strat run with F r = 0.168, with the lowest values of N and f . For other runs shown, while minor deviations are evident, the spectra by large remain close to k −5/3 scaling.
On the other hand, in Fig. 2b, the spectra for runs with stronger stratification, corresponding to F r ≤ 0.069 (or equivalently R IB < 1) are shown. The spectrum for the HIT run (dashed black line) is also included for comparison. As evident, this regime significantly differs from that at weaker stratification, with the spectra exhibiting a steeper spectral slope than k −5/3 . It has been observed in the closely related context of a time-dependent flow [28,42] that the spectra followed the Bolgiano-Obhukov prediction [43,44]. In this regime, energy transfer is strongly affected by the temperature fluctuations (potential energy). The straight line in Fig. 2b indicates the corresponding k −11/5 prediction. Given the small range of spatial scales available in our simulations, however, it is difficult to conclusively demonstrate that the spectra shown indeed conform to a −11/5 scaling. Nevertheless, the sharp change in spectral slopes clearly evidences that the energy transfer between scales is strongly perturbed in runs with R IB < 1, and much less so in runs with R IB > 1. In fact, this demarcation based on R IB will become increasingly prominent as more results are discussed in later sections. The spectra of the temperature fluctuations, although not shown, follow similar trends as that of velocity. We note that the interesting possibility of having a constant-flux inverse cascade of energy towards large scales, together with a constant-flux direct cascade of energy towards small scales, induced by the combined effect of RaS [24,45], does not seem applicable to the runs studied here. Indications of inverse cascade, as evidenced by an energy flux towards scales larger than that of the forcing, were reported in flows with comparable Rossby and Froude numbers, albeit with larger buoyancy Reynolds numbers than here [24,34,45]. In our runs, the lack of energy accumulation at scales larger than L (which equals the forcing scale L f ) is likely attributable to the lack of scale separation between L and the length scale of the flow domain.
The above transition between strong hydrodynamic turbulence and regimes which are strongly dominated by RaS, as evidenced from the energy spectra, clearly warrant further inspection. An important aspect of the flows in this study is the existence of a strong anisotropy of velocity fluctuations. Namely, the presence of RaS induces very different properties of the velocity fluctuations in the horizontal plane and in the vertical direction. In this regard, we consider the components of the velocity u ⊥ and u , respectively in the directions perpendicular and parallel to gravity. The u ⊥ is obtained by averaging over the two components in the x-y horizontal plane (which is also the plane of rotation) and u simply corresponds to the vertical z-component. An obvious manifestation of the anisotropy between the two directions is in variances of u ⊥ and u , as shown in Table II. It is evident that apart from the HIT run (case 0) and to large extent the rotation only run (case 8), the variances for both components are very different. To precisely quantify the anisotropy, we consider the ratio ξ u = [ u 2 ⊥ / u 2 ] 1/2 , also shown in Table II. Expectedly, this ratio is unity for case 0 and very close to unity for case 8. However, for all runs with stratification, the ratio is always greater than 1. For cases 1-3, the ratio increases slowly with decreasing F r, whereas a sharp increase occurs for cases 4-6 (and also case 7). Clearly, this transition corresponds to the change in spectral slope observed for the spectra in Fig. 2. The runs with smaller ξ u correspond to spectra which show k −5/3 scaling, whereas runs with sharp increase of ξ u correspond to spectra with steeper Bolgiano-like scaling. In these flows, the fluctuations of u are essentially suppressed when imposed RaS are strong. We note that the values shown here are consistent with some earlier reported runs in [25]. Variance and kurtosis of velocity components. We use ξ u = u 2 ⊥ 1/2 / u 2 1/2 , to measure the anisotropy, whereas K u and K u ⊥ are the kurtosis of u and u ⊥ respectively. The skewness of these components (not shown) is consistent with 0. Given the difficulties in determining the kurtosis K u for run 6 (F r = 0.045), we stress that the accuracy of the indicated value -marked with an asterisk -is low.
To further investigate the anisotropy, we next consider higher order moments of the velocity components. Since the probability density functions (PDFs) of individual velocity components are symmetric (due to underlying symmetry of Boussinesq equations), we note that the values of the third moments of the velocity fluctuations are expected to be zero. Our results, although not shown here, are compatible with this. Thus, we consider the fourth-order moment, through the kurtosis, defined as It is known that the PDFs of individual velocity components in HIT are approximately Gaussian, which implies that the kurtosis is approximately 3. Consistent with that, we find that the value of the kurtosis is very close to 3 in the case HIT, (see Table II). In the case of a purely stratified flow, however, it was found that the distribution of the velocity component parallel to the direction of stratification, u , could be significantly wider than Gaussian [32,33]. The strong deviations from a Gaussian distribution have been shown to originate from intermittent vertical drafts [33]. On the other hand, the kurtosis of u ⊥ is always found to be approximately 3. Our values, as shown in Table II, are consistent with these trends. We find that K u ⊥ is slightly smaller than 3 at low stratification, and slightly larger than 3 for runs with strong stratification (F r 0.08), with PDFs (not shown) being approximately Gaussian. In the vertical direction, K u increases very significantly beyond 3 as strength of RaS increases.
It is interesting to contrast the dependence of K u in our runs to the stratification-only runs reported in [33] (see in particular their Fig.5). While the sharp transition in K u appears to be qualitatively similar. However, quantitatively the large values of K u for RaS runs appear to occur at lower F r values than in [33]. For even smaller F r (not considered in this work), in wavedominated regime, one can expect that K u for RaS runs will again become 3 [33], as it was the case for the stratification-only run (case 7) discussed here. A more systematic quantification of K u for RaS runs will be addressed separately.

IV. ACCELERATION STATISTICS
The acceleration experienced by a fluid element, defined by the rate of change of velocity in the Lagrangian frame, i.e., a + = du + /dt and resulting from the balance of forces acting on it, is arguably the simplest descriptor of its motion, as also directly reflected in the governing fluid equations. In addition to its fundamental importance in turbulence theory [12], a key motivation to study acceleration comes from its central role in stochastic modeling of turbulent dispersion [46]. While acceleration statistics have been thoroughly investigated in isotropic turbulence, they still remain poorly understood in rotating-stratified flows. In the following sub-sections, we study different aspects related to acceleration, namely, acceleration variance and kurtosis, the probability density functions and Lagrangian autocorrelations and frequency spectra. A central theme is to  Fig. 4). The underlying symmetry of the problem imposes that the skewness (not shown) is zero. The values of kurtosis marked with an asterisk correspond to relatively low accuracy.
understand the effect of imposed anisotropy. Hence, similar to our earlier analysis of velocity field (in Sec.3), we study the properties of Lagrangian acceleration by considering separately the horizontal and the vertical components, i.e., a + ⊥ and a + respectively. For convenience, we will omit the superscript '+' henceforth, since our analyses will only involve Lagrangian quantities (unless otherwise specified).

A. Acceleration variance
Extending the earlier analysis of the Eulerian velocity field, we first examine the anisotropy in acceleration by considering the second and fourth order moments, namely, the variance and kurtosis respectively. The moments are listed in Table III. Due to the underlying symmetry of the Boussinesq equations, the third and all other odd moments of acceleration components are zero (similar to the velocity components). As can been seen from cases 1-6, the properties of a ⊥ and a show starking differences. While the variances of a rapidly decrease with decreasing Ro (and F r), the variances of a ⊥ do not change that much. To better understand this behavior, we consider the anisotropy ratio ξ a = [ a 2 ⊥ / a 2 ] 1/2 . For HIT run (case 0), the ratio, as expected, is equal to 1, whereas for rot-strat cases two separate patterns are visible.
The first corresponds to cases 1-3 with relatively weak rotation and stratification (RaS), where the ratio ξ a is approximately constant, and slightly smaller than 1. The ratio being close to unity suggests that the effect of anisotropy on small scales is only minor. Furthermore, a tentative explanation for the ratio slightly smaller than unity can be proposed, based on the observation that the relative strength of stratification is significantly stronger compared to that of rotation (given N/f = 5). For that reason, the acceleration variance in the vertical direction is likely to experience more enhancement due to stratification than that in the horizontal direction.
On the other hand, runs 4-6 correspond to significantly stronger RaS, which leads to a different pattern, with ξ a > 1 and further increasing with decreasing F r (and Ro). Since runs 4-6 correspond to R IB < 1, they are in a regime where waves play a stronger role, compared to (turbulent) eddies. With strong stratification, the motion in the vertical direction is strongly suppressed, as also reflected in ξ u values (see Table II). With relative weaker effect of rotation, the acceleration in the horizontal direction is still affected by turbulent eddies, resulting in ξ a to become greater than 1. Thereafter, the effect becomes even more pronounced with further increase in strength of RaS.
In this regard, further qualitative insight can be obtained by considering the ratio ξ u /ξ a , also shown in Table III. For HIT, this ratio is unity, but with increasing strength of RaS, we observe ξ u /ξ a increases, peaking at around 6 for case 4, and thereafter decreasing slightly with further increase in strength of RaS. A limiting value can be obtained in the case where RaS is very strong, such that the flow is dominated by linear processes (with turbulent eddies playing a very weak role). This would lead to ξ u /ξ a ≈ N/f , by assuming that a ⊥ ∼ u ⊥ f and a ∼ u N . The value for case 6, with ξ u /ξ a = 18.4/3.8 ≈ 5, appears to be consistent with this consideration, though it remains to be further tested for different N/f values.
It is instructive to compare the acceleration variances for case 5 (with strong RaS), to case 7, with exactly the same stratification, N , but no rotation and to case 8, with the same rotation, f , but no stratification. The variances of a ⊥ and a are stronger in the presence of both RaS, than with either of the two effects taken separately. The anisotropy ratio ξ a also shows a very striking difference, when comparing runs 5, 7 and 8. For rotation-only (run 8), we observe ξ a ≈ 1, implying little to no small-scale anisotropy. This can be viewed as consistent with the observation that the flow spectrum for run 8 is overall comparable to the Kolmogorov spectrum (see Fig. 2a), so the flow is dominated by eddies. On the other hand, the anisotropy ratio is ξ a < 1 (0.66) for run 7, implying stronger acceleration variances in the vertical than in the horizontal direction, but an opposite situation in the presence of both RaS: ξ a > 1 (2.1). The difference can be traced back to the much larger value of a 2 ⊥ in the presence of rotation, a consequence of the enhancement of the horizontal component of acceleration due to the Coriolis force.
It is interesting to notice, when comparing the anisotropy ratios ξ u and ξ a for the velocity and acceleration components, that ξ a is generally much closer to unity than ξ u . This is consistent with the expectation that acceleration is sensitive mostly to the small-scales of turbulence, which are, on general grounds, close to locally isotropic, the more so as the Reynolds number Re increases [41]. To further explore this aspect, we consider the following relation from classical turbulence which results from the application of Kolmogorov's similarity hypothesis to acceleration variance [7]. Here, a 0 is a dimensionless constant, which is plausibly universal for HIT. Empirical evidence based on studies of isotropic turbulence has shown that a 0 increases slowly with Reynolds numbers [47], which can be possibly viewed as a manifestation of small-scale intermittency, unaccounted for in Kolmogorov's 1941 theory. In contrast, for rotating-stratified turbulence, one may expect additional deviations because of imposed anisotropy, especially if the Re is not very large. Fig. 3 shows the values of the a ⊥ 0 and a 0 as a function of F r for cases 1 through 6, together with the value obtained for case 0, in the absence of RaS. The values of a 0 and a ⊥ 0 do not vary very much for 1/F r 12.5, and sharply increase, especially in the perpendicular direction, at large values of 1/F r. In particular, the value of a ⊥ 0 at 1/F r ≈ 16 exceeds that in the HIT case by a factor larger than ∼ 10. This points to a strong difference of the small-scale properties between runs 4-6 and the HIT run (case 0). In comparison, and despite the larger value of N (N/f = 5 for runs 1-6), the value of a 0 remains comparable to that obtained in the HIT case. In the purely stratified case (run 7), the value of a 0 is slightly reduced compared to case 5 in the presence of both RaS (a 0 ≈ 3.1). On the other hand, the value of a ⊥ 0 ≈ 1.4 for case 7 is reduced by over an order of magnitude compared to case 5. For the rotation-only run, a 0 is not very different compared to the HIT case, suggesting that the small scales are still dominated by turbulence, even though the large scales are highly anisotropic. This demonstrates that acceleration becomes progressively dominated, as Ro and F r decrease, by the frequencies associated with rotation/buoyancy. As a consequence, the values of a 0 shown in Fig. 3 do not reflect only the small-scale structure of the flow, but also its global properties.
It is useful here to notice that the various dimensionless numbers characterizing the flow, defined by Eq. (5) provide a quantitative measure of the physical effects influencing the flow. In particular, we recall that the values of R IB directly compare the mean dissipation rate and the Brunt-Väisälä frequency. As shown in Table I, R IB is larger than 1 for runs 1-3, which indicates that the vorticity fluctuations prevail over N . In contrast, R IB is smaller than 1 for runs 4-6, pointing to a dominance of the waves. This distinction coincides with the notion that the acceleration variance is dominated by vorticity, in runs 1-3, whereas it is rather due to waves in runs 4-6. The other (buoyancy) Reynolds numbers, R B , and the micro-Rossby number, Ro ω , provide a consistent picture, although the values of the numbers are not as simple to interpret.

B. Acceleration PDFs
The analysis of acceleration variance can be generalized by considering the probability density functions (PDFs) of acceleration components. In HIT, similar to Eulerian velocity gradients, acceleration is also characterized by extremely large fluctuations reflected in broad tails of the PDF, which further broaden with increasing Reynolds number [12,35]. In this subsection, we study the effect of imposed RaS on the PDFs of acceleration, especially focusing on the long tails. Once again, we separate out the contributions a and a ⊥ . Fig. 4 shows the standardized PDFs of a ⊥ (left column) and a (right column) for all the runs listed in Table III. The cases with R IB > 1 are shown in the top row, whereas cases with R IB < 1 are shown in the bottom row. Two very distinct behaviors can be observed once again. At moderate RaS (top panels) the PDFs of acceleration exhibit very broad tails, and in this sense, they differ only quantitatively from those obtained in the HIT case. This property is also reflected by the values of the kurtosis of the distributions, shown in Table III, which are all 20, even larger than the values obtained for the HIT case. The run corresponding to pure rotation also exhibits very high kurtosis of the acceleration. Note that the accuracy of the estimates of kurtosis indicated in the table are of the order of 10% or less, except for values marked with an asterisk, which exhibited very large fluctuations, leading to higher error bars (of at least 20 − 30%).
In contrast, the bottom panels of Fig. 4 shows that at stronger stratification the tails of the PDFs become significantly less extended. As shown in Table III, the corresponding kurtosis factors are down to values smaller than in the HIT run, although still larger than Gaussian value of 3. This is particularly clear for the horizontal component a ⊥ of acceleration, see Fig. 4c, and to a lesser extent, for the vertical component a (Fig. 4d). The tendency of the PDFs to become narrower is common to the runs with a strong stratification; rotation affects only quantitatively the PDF. The observation that, when the stratification is strong, the large values of acceleration are reduced (see Fig. 4), is generally consistent with the values of the kurtosis K a and K a ⊥ , listed in Table III. It is worth recalling that for runs 4-6, the values of a 0 were found to be much larger than what is expected, based on the eddy motion only (see Fig. 3), therefore pointing to a strong role of waves in the motion of tracers. We view the two facts as a manifestation of the strong role of waves for runs 4-6. In addition, we notice that the PDF of a , normalized by its rms are consistently broader than those of a ⊥ . In a flow dominated by waves, this qualitative difference may appear as surprising. Even in the regimes with strong RaS, however, nontrivial hydrodynamic effects play an important role [33]. We attribute the comparatively stronger values of K a to the strong bursts of vertical velocity, clearly documented in purely stratified flows [33]. These bursts leave a strong signature on  Table II. The intermittent episodes of strong vertical motion are leading to rare, intense vertical acceleration events, clearly seen in Fig. 4.

C. Acceleration autocorrelation and frequency spectra
Essential information on the motion of tracers can be obtained from the autocorrelation function of quantities fluctuating along Lagrangian trajectories. For example, the autocorrelation function can be used to determine relevant time scales and to form the frequency spectrum via Fourier transform. Since our analysis is concerned with statistically stationary signals, the autocorrelation only depends on the chosen time lag. The autocorrelation ρ a (τ ) is defined as C a (τ ) = a(t + τ )a(t) , ρ a (τ ) = C a (τ )/C a (0) (9) where τ is the time lag and C a is the autocovariance. Statistical stationarity implies ρ a (−τ ) = ρ a (τ ) and also imposes that the integral of the acceleration autocorrelation must be zero, i.e., ∞ 0 ρ a (τ )dτ = 0 [48]. The frequency spectrum of acceleration E a , which is the Fourier transform of C a (τ ) can be defined as By definition, the integral of the spectrum gives the variance of the signal. Yet again, we consider a ⊥ and a separately. Perpendicular component: Fig. 5 shows the autororrelation function ρ a ⊥ (τ ) for all the runs with RaS. In Fig. 5a, the time lag is normalized with τ η , the Kolmogorov time scale characterizing the small-scale motion. We also include the HIT case (dashed line) and the pure rotation case (dash-dotted line) as a reference. For HIT, the acceleration autocorrelation is known to decay rapidly, becoming zero at τ ≈ 2τ η and thereafter becoming negative and slowly approaching zero again at very large time lags [49]. This is clearly seen in Fig. 5a, see the dashed line corresponding to our HIT run. Interestingly, we observe that the runs with relatively weak RaS, i.e., F r ≥ 0.114,  and also the run with only rotation, behave very similarly to the HIT case (see inset of Fig. 5a), albeit with minor variation in the zero-crossing point. For the rot-strat run with F r = 0.086, the autocorrelation function shows a similar behavior up to very small time lags, but thereafter deviates with the zero-crossing extending to τ ≈ 8τ η . The large time behavior is markedly different with very low frequency oscillations. On the other hand, for runs with strong RaS, i.e., F r ≤ 0.07, a very different behavior is observed right from small time lags. The autocorrelation shows strong oscillations which are eventually damped out at large time lags. Furthermore, the zero-crossing point also is strongly dependent on F r, moving to smaller time lags with decreasing F r. Plotting ρ a ⊥ as a function of τ f /2π, see Fig. 5b, shows that the period of the damped oscillations is close to 2π/f (with small but significant deviation). For comparison, we also show in Fig. 5b the function ρ a ⊥ for the case of a purely rotating flow. Somewhat surprisingly, the tendency to oscillate at a frequency ∼ f is not as strong as in the case where both RaS are present. A similar trend was observed in [13]. We notice, however, that the stronger stratification (N = 5f in runs 1-6) implies the presence of faster waves than in a purely rotating flow at the same Ro number, involving in particular high values of k ⊥ , since the dispersion law for inertia-gravity waves can be written as k 2 ω 2 IG = N 2 k 2 ⊥ + f 2 k . To better identify the frequencies present in the autocorrelation function, it is apt to take its Fourier transform to obtain the frequency spectrum. Fig. 6 shows the spectrum E a ⊥ . The spectra have been made dimensionless by dividing by the mean dissipation rate ǫ . In Fig. 6a, the frequency is normalized by τ η , whereas in Fig. 6b, we use the rotation frequency f . Note that the runs shown in Fig. 6a and b correspond to the same runs shown respectively in Fig. 5a and b. For the HIT run in Fig. 6a, Kolmogorov scaling would imply that the spectrum E a is constant in the inertial range (defined as τ η /T E ωτ η 1). While the validity of Kolmogorov's similarity hypotheses to Lagrangian statistics is debatable, DNS studies suggest that the above scaling law, i.e., E a ∼ ǫ , still may be approximately satisfied [50]. Consistent with this, we find that the spectrum E a ⊥ (ω) is approximately flat over a limited range of values of ω. With addition of weak RaS, the changes in the structure of the spectrum appear to be relatively minor. In particular, for runs with F r ≥ 0.086, the differences in spectra, relative to the HIT case, only become significant for ωτ η 0.2. For the particular case of F r = 0.086, a small peak emerges around ωτ η ≈ 0.1. With increasing strength of RaS, the differences become more pronounced, with all the cases for F r ≤ 0.086 now showing a prominent peak (which is relatively broad). This observation is consistent with the behavior of the autocorrelation function at small times, as observed in Fig. 5. In addition, the increase in values of E a ⊥ / ǫ also corresponds to the strong increase in the value of a 0 for small F r cases, as shown in Fig. 3 (for the perp. data points). Interestingly, the decay of spectra for all cases are somewhat similar at very large frequencies (and not very different from the HIT case). This possibly confirms that the role of turbulent eddies in the horizontal direction is still relevant as compared to vertical direction, at least in the parameter range covered in this work (note once again the fact that N/f is 5, plays an important role, since it renders the effect of rotation relatively weaker than that of stratification) In Fig. 6b the spectra are shown as a function of ω/f . As expected, the peaks observed at low frequencies in Fig. 6a, are all centered around ω/f ≈ 1. For the runs with strong RaS, the frequency corresponding to stratification, i.e., ω/f = N/f = 5, does not seem to play any particularly important role.
Parallel component: While the role of rotation was dominant in the perpendicular direction, we can in turn expect stratification to dominate in the parallel direction. In Fig. 7a and b, we show the autocorrelation in the parallel direction ρ a as a function of τ /τ η and of τ (N/2π) respectively. As it was the case for the perpendicular component of acceleration, the deviation of the autocorrelation function from the reference homogeneous isotropic case, shown as a dashed line, increases gradually when F r decreases. However, contrary to ρ a ⊥ , ρ a tends to decay significantly faster than in the HIT case when F r diminishes. In addition, for all cases, overshooting oscillations are clearly visible, with the amplitude particularly strong for cases with F r ≤ 0.086. While the zero-crossing of ρ a for cases with weak RaS still appears to be close to τ ≈ 2τ η (with minor deviations), for cases with strong RaS, the zero-crossing shifts to even smaller time lags of τ ≈ 0.5τ η . This is in contrast to ρ a ⊥ , where the zero-crossing first shifted to higher time lags, and then moved to smaller time lags with decreasing F r. Fig. 7b clearly reveals that the period of the oscillations is equal to the Brunt-Väisälä period, 2π/N . This strongly suggests that the motion in the parallel direction is dominated by stratification, and that turbulence plays a much smaller role for the vertical motion (in comparison to horizontal motion, where even at strongest rotation rate, the role of turbulence still could not be ignored). For comparison, the autocorrelation function ρ a is shown in Fig. 7b for the run with a strong stratification, F r = 0.03. The tendency to oscillate is even stronger in this run, an effect amplified by the low value of F r (and R IB ) in this case (see Table II). The corresponding frequency spectra E a (ω) are plotted as a function of ωτ η and of ω/N in Fig. 8a and b respectively. The spectra are again non-dimensionalized by ǫ . Unlike in the case of E a ⊥ , we observe that deviations from the HIT run are already prominent even for runs with weak RaS. All spectra are characterized by the presence of a peak, with the peak becoming more sharp and prominent as F r decreases. An inspection of spectra in Fig. 8b clearly shows that all these peaks correspond to ω/N ≈ 1, i.e., the respective Brunt-Väisälä frequencies, conforming with the dominant role of stratification in vertical direction. Note the peak for stratification-only run is even more pronounced. For runs with strong RaS, i.e., F r ≤ 0.069, a very minor enhancement of the spectra is visible in the band of frequencies ω/N ≈ 0.2 (or ω/f ≈ 1).
Another point to note is that even though the peaks in E a become sharper with decreasing F r, they are all still at comparable values, in contrast to E a ⊥ , where the peaks also increased in value by more than an order of magnitude. This can be explained by the variation of acceleration variance as seen in Fig. 3. While the a ⊥ 0 sharply shoots up for runs with F r ≤ 0.069, the a 0 only shows a minor variation in comparison.

D. Discussion
The analysis presented in this section suggests that the behavior of acceleration statistics, studied here for N/f = 5, and over a limited range of Reynolds number, can be separated into two qualitatively different regimes. The first one corresponding to weak RaS, such that the rms of the vorticity fluctuations is larger than the Brunt-Väisälä frequency: R IB > 1. In this regime, the acceleration properties are close to those documented in HIT. The clear large scale anisotropy of the flow, as seen from the velocity statistics listed in Table II, only moderately affects the acceleration properties. This is consistent with the observation, see Fig. 2a, that the spectrum in this regime is close to that predicted by Kolmogorov theory.
In the second regime, where the vorticity fluctuations are smaller than the Brunt-Väisälä frequency (R IB < 1), the statistical properties of accelerations are totally different, as strong RaS clearly affect the small-scales. This is in full qualitative agreement with the flow spectra, see Fig. 2b, which deviate very strongly from those predicted by Kolmogorov theory. Although one may have surmised that in this regime where waves dominate, the acceleration statistics are close to Gaussian, we observe, surprisingly, that the kurtosis of acceleration components are much larger than 3 (with a stronger kurtosis for the vertical, than for the horizontal component of acceleration). This unexpected phenomenon can be related to the observation of bursts of vertical velocity, manifested by the large values of K u listed in Table II, particularly for runs 4-6. These bursts are related to the phenomenon observed in [33] in the case of purely stratified flows, in comparable domains of parameters. We expect that the mechanisms are similar in presence of a relatively weak rotation (N/f = 5), as it is the case in the present study.
We stress that the picture we are drawing is based on runs at a relative moderate Reynolds numbers. A natural question arises about whether these conclusions would apply at even larger Reynolds numbers. Further investigations, varying not only the Reynolds number, but also the ratio N/f , should bring more insight on the structure of these flows.
Whereas R IB seems, in the present problem, to provide a clear criterion to distinguish between the different regimes observed here, alternative ways [25] have been proposed to describe the transition between wave-dominated flows, and flows where waves and eddies interact. Namely, it was found that the transition occurs around R B ∼ 10. In this sense, the flows studied here are all in the regime of interacting waves and eddies, although run 6 is very close to the transition towards a wave-dominated flow.

V. LAGRANGIAN VELOCITY AND DISPERSION STATISTICS
The spread or dispersion of a material under the action of turbulence is of obvious importance and motivates the study of Lagrangian velocity and dispersion statistics. In this regard, we briefly summarize the classical theory below [7].
Single-particle dispersion is best understood by considering the mean-square displacement of a particle from its initial position. From Eq. (7), we can write is the displacement of the particle from its position x i (0) at t = 0 to its position x i (t) at time t, and ρ u i (t ′ , t ′′ ) is the velocity autocorrelation. For statistical stationarity, the autocorrelation function at t ′ and t ′′ only depends on τ = |t ′ − t ′′ |, allowing Eq. (11) to be rewritten as Once again, we treat separately the horizontal and vertical components, as explained before. The expression for the dispersion, Eq. (12), simplifies in the limit of short and long times, respectively to where D ⊥, = 2 u 2 ⊥, T ⊥, L are the diffusion coefficients and (15) are the Lagrangian integral times. In the following subsections, we discuss separately the behavior of the autocorrelation functions, integral time scales and then the mean-square displacements.

A. Velocity autocorrelation and spectrum
In this subsection, we investigate the behavior of the velocity autocorrelations and the corresponding frequency spectra. We recall that the velocity frequency spectrum can be simply derived from the acceleration spectrum using the relation E u (ω) = E a (ω)/ω 2 [48]. Thus in that regards, E u does not provide completely novel information. However, in flows considered in this work, where the role of the frequencies associated with rotation and stratification (RaS) needs to be quantified, looking at both E u and E a is useful to understand how velocity and acceleration are affected.
Perpendicular component: Fig. 9a-b shows the autocorrelation functions ρ u ⊥ (τ ) along with the corresponding frequency spectra E u ⊥ (ω) in Fig. 9c-d. Similar to results for perpendicular component of acceleration, the time lag and frequency are normalized first by τ η and then by f . The correlation function ρ u ⊥ , plotted as a function of τ /τ η , see Fig. 9a, shows only weak deviations from the HIT case (shown in black dashed line) when RaS are moderate (runs 1-3). We recall that in the HIT case, the long time behavior of the velocity autocorrelation can be well represented by an exponential functional form, ρ u ∼ exp −τ /T L [11]. However, the deviations become stronger with imposed RaS. In particular, the runs with RaS are fundamentally different with appearance of oscillations resulting in negative values of the correlation function.
To better understand these oscillations, Fig. 9b shows the autocorrelation as a function of τ f /2π. Although an oscillatory behavior can be seen in the runs with the highest values of RaS rates, the frequency is close to, but differs from f . On the other hand, ρ u ⊥ in the flow with rotation only (shown in black dash-dotted line in Fig. 9b), appears to be monotonically decreasing (no oscillations). Further insight can be obtained from the frequency spectra, shown in Fig. 9c-d. Fig. 9d shows the spectra as a function of ω/f , and reveals, in the presence of RaS, the presence of a weak maximum, in a band of frequencies around ω/f ≈ 1. This differs from the acceleration autocorrelation function, where the rotation frequency was more prominent with increasing strength of RaS. Fig. 9c shows the spectra as a function of ωτ η , with two separate behaviors clearly evident. At weak RaS (F r 0.08), the spectra are very similar to the HIT run, with only small deviations at low frequencies. At stronger RaS (F r 0.07), the spectra increasingly differ at low and moderate frequencies, with the emergence of strong activity in a broadband around ω/f = 1 as identified in Fig. 9d. At high frequencies, on the other hand, all the spectra approximately collapse and are close to HIT. These observations, similar to those for acceleration frequency spectra, once again suggest that the role of turbulence in the horizontal plane is always significant, at least over the parameter range considered in this work. This is consistent with the observation that Ro ω > 1 for all our runs, suggesting that turbulent eddies always prevail over rotation.
It is interesting to note that in the case of pure rotation, the frequency ω = f does not seem to play any particular role in the spectrum, even less than in the presence of stratification. We also notice that the value of Ro ω is actually lower for run 5 (in the presence of stratification) than for run 8 (without stratification), both at f ≈ 2.95. These observations indicate that the  Fig. 10. We first consider the autocorrelation ρ u as a function of τ /τ η in Fig. 10a. In all the cases with RaS, ρ u strongly differs from the HIT case (indicated by a black dashed line). Two patterns are distinctly visible. First, all the autocorrelation functions decay rapidly. In terms of the dimensionless time t/τ η , the decay becomes increasingly rapid when F r decreases. Second, the autocorrelations all overshoot to become strongly negative and show distinct oscillations. The primary period of these oscillations corresponds to the stratification frequency, as demonstrated by Fig. 10b, which shows the autocorrelations as a function of τ N/2π. In addition, Fig. 10b also shows that all cases superpose reasonably well for the initial decay of the autocorrelation, These results clearly demonstrate that stratification plays a dominant role in the vertical motion, and in comparison turbulence plays a minor role. This is true, even for the runs at weak RaS (F r 0.08). A similar behavior was obtained for the vertical acceleration correlation function, see Section IV C. We notice that in the case of a purely stratified run, without any rotation, the autocorrelation function ρ u clearly oscillates, with a modulation at twice the Brunt-Väisälä period.
To complement our investigation of the autocorrelation, Fig. 10 shows the frequency spectra E u as a function of ωτ η (panel c) and ω/N (panel d). The deviation from the spectrum E u for the homogeneous, isotropic flow, shown with dashed lines in panel c, clearly increases when Ro diminishes in the range ωτ η ≤ 1. The increase in the Brunt-Vaäisälä frequency leads to the formation of a relatively narrow peak. Panel d demonstrates that the peak is centered at N , with a weak modulation at the lower (Coriolis) frequency, f = N/5. In the case of a flow with stratification only, no rotation (indicated by a dashed-dotted line in panel d), the sharp peak at ω ≈ N is accompanied by a secondary peak at frequency ∼ N/2, consistent with the observation of ρ u in Fig. 10b

B. Integral time scales
The Lagrangian integral time scale T L , which characterizes particle dispersion at large times, is simply obtained by integrating the correlation function: T L = ∞ 0 ρ u dτ . For this reason, we discuss first T L before we present the results on mean-square displacement. In addition to T ⊥ L and T L , we also introduce T all L , defined based on all three components of velocity, which characterizes the total dispersion. We begin by noticing that many of the correlation functions shown in Figs. 9 and 10 exhibit an oscillatory behavior. This makes the calculation of their integrals generally prone to statistical errors. As such, the values given in this section are only indicative. The integral time scales for all the cases are listed in Table IV. For the horizontal (perpendicular) direction, the integral time scale is compared with τ η and also with the time period 2π/f associated with rotation. For the vertical (parallel) direction, the time period 2π/N is used instead.
The variation from run to run of the time T ⊥ L in the horizontal direction, is found to be rather weak, at least at moderate stratification, rotation and Reynolds number. In particular, the ratio T ⊥ L /τ η , is of the order of 10 for the HIT run, and remains approximately constant for cases 1-4 and also for the run 8, which corresponds to a flow with rotation only. For cases 5-6, the ratio drops by a factor ∼ 2, which can be attributed to a sharp decrease in the mean dissipation ǫ , leading to an increase in τ η . The above observations are consistent with earlier results which also highlight a sharp transition to wave dominated regime as R IB significantly decreases. We note that for the run with pure stratification (case 7), the velocity autocorrelation in the horizontal direction converges extremely slowly to zero at large time lags, pointing to an extremely long value of T ⊥ L . A similar behavior has been noticed in [19]. In this regard, the presence of rotation plays a crucial role. In particular, the dispersion relation and the structure of the eigenmodes, underlying the wave motion is strongly modified by rotation. On the other hand, we find that the ratio T ⊥ L f /2π increases as the rotation and stratification increase.
We find the integral time in the vertical direction, T L , to be substantially smaller than T ⊥ L for all runs with stratification (cases 1-7). The ratio T L /τ η is still around 10 for the HIT and purely rotating flows (cases 0 and 8 respectively), but is greatly reduced for all runs with stratification. This can be evidently attributed to the strongly oscillating nature of the velocity autocorrelation functions for cases 1-7 (see Fig. 10b), which results in significant cancellation when calculating T L . The ratio T L /τ η drops by more than an order of magnitude in going from HIT to case 1, which corresponds to the weakest RaS. Thereafter, with increasing strength of RaS, the ratio further decreases. Interestingly, the ratio T L N/2π shows a much reduced variation, and remains in the range 0.01-0.02 for all cases with stratification. This possibly suggests that despite the strong oscillations in the autocorrelation the time scale, T L , is not strictly zero, but rather scales as N −1 , with a very small numerical prefactor. Finally, as a consequence of T L ≪ T ⊥ L we find T all L ≈ T ⊥ L . As the ratio T all L /τ η does not vary very much, plotting the velocity autocorrelation, ρ u ⊥ , as a function of τ /T ⊥ L leads to a Figure qualitatively very similar to Fig. 9a.
C. Mean-square displacement Fig. 11 shows the mean-square displacement for horizontal and vertical directions as a function of time (both axes normalized appropriately by Kolmogorov scales). At small times, a clear t 2 scaling is visible for for all runs in both components, as expected from Eq. (13). The displacement in the horizontal direction, Y 2 ⊥ (t), is enhanced in the vertical direction by RaS, whereas it is reduced in the vertical direction. This can be readily explained by considering the variance of the components of the velocity v 2 ⊥ and v 2 , as shown in Table II (see also Eq. (13)). At long times, transition to a diffusive regime, given by Eq. (14), is expectedly seen in the HIT run [48], for both horizontal and vertical components. A similar behavior is also visible for runs with weak RaS (F r 0.08), when considering the horizontal component. Whereas the growth of the vertical component for these runs is significantly inhibited in comparison, although a slow approach to the diffusive regime is nevertheless noticeable. Note the rotation-only run, behaves like the HIT run, consistent with the small degree of anisotropy in both velocity and acceleration statistics as noted earlier.
On the other hand, for runs with strong RaS (F r 0.07), the growth of both horizontal and vertical components appears to be slower than a linear behavior. This is particularly clear for the vertical displacements Y 2 , which appear to increase extremely slowly with time, if at all. In fact, the values of Y 2 /η 2 remain small, of the order of ∼ 10, over a very long time. The absence of a clear diffusive regime for Y 2 is arguably not very surprising, given the small values of T L and v 2 , hence of the diffusion coefficient D = v 2 T L . Assuming that T L behaves as 2π/N (multiplied with a small prefactor), as suggested by the results shown in Table IV, the expected diffusive behavior will be visible only after an extremely long time (possibly a factor 10-100 longer than used in present work). Consistent with our own results, the near constancy of Y 2 was also reported in purely stratified flows [19,51].
Contrary to the very low values of T L for runs with strong RaS (see Table IV), the values of T ⊥ L do not appear to be particularly small. This makes the lack of a clear diffusive regime for runs 4-6 for the horizontal displacement much more surprising. It has been noticed in the context of stratified flows, that the horizontal displacement could grow as t 2 [51]. Here, our results show that, in a regime dominated by waves (runs 4-6), the growth of Y 2 ⊥ is slower than t 1 over a long time. This suggests a strong interaction between the horizontal and vertical directions owing to the presence of rotation (and inertia-gravity waves), which was also visible in autocorrelations of acceleration and velocity. However, given the large values of T ⊥ L , it can be expected that the diffusive regime will ultimately be reached at sufficiently longer times.
Finally, we note that the strong vertical drafts responsible for high kurtosis of u (see Section III), and increased mixing efficiency [33], do not appear to cause an appreciable growth in Y 2 for the given set of runs considered here. This is not inconsistent given the kurtosis is a normalized fourth-order moment, whereas Y 2 goes as the variance of vertical velocity and the integral scale T L , both of which are strongly suppressed. It is likely that at higher Reynolds numbers than considered here, the affect of the vertical drafts will be felt more directly on the vertical displacement of fluid particles.

VI. DISCUSSION AND CONCLUSIONS
Motivated by the well established effect of both rotation and stratification (RaS) in numerous geophysical flows, and by the observation of Lagrangian tracers, such as buoys in the ocean [52], we have investigated here the properties of particle trajectories in such flows. We have utilized direct numerical simulations of Boussinesq equations in a periodic domain, with both stratification and the axis of rotation aligned in the vertical direction. We have focused on cases where the corresponding rotation and stratification frequencies, f and N respectively, are held in a fixed ratio N/f = 5, as relevant in oceanographic situations [21]. The usual computer limitations have restricted the Reynolds numbers to approximately around 4000 in our runs, which are smaller than those observed in nature, but still large enough to allow turbulence to adequately develop. Whereas the Froude numbers investigated are in the range 0.03 F r 0.2, with corresponding Rossby number Ro = 5F r, compatible with geophysical fluid flows. Our results illustrate the complex physical effects involved in a turbulent flow under the combined action of RaS. Whereas the corresponding physical effects have been recently investigated in an Eulerian context in many studies, we have documented here the elementary properties of particle trajectories, in particular concerning acceleration, velocity and displacement statistics. Due to the imposed large-scale anisotropy, we differentiate between motion in the horizontal plane and in the vertical direction. We also appropriately compare with rotation-only and stratification-only runs and also with the relatively well studied case of isotropic turbulence.
While the parameters F r and Ro vary smoothly with the strength of imposed RaS, we observe a sharp change in the nature of the flow when F r 0.07 for the runs considered, corresponding to a regime where waves appear to dominate over the non-linearities, as also observed in [33] for stratification-only runs. An inspection of the Eulerian energy spectra (and the temperature spectra), shows a clear transition from a Kolmogorov k −5/3 like behavior to a Bolgiano k −11/5 like behavior at intermediate wavenumbers. The transition also corresponds to the peak of intermittent disruption of shear layers, with a maximum of the kurtosis of the vertical velocity. In this context, our work extends the stratification-only results of [33]. Interestingly, we find that this transition to wave-dominated regime can be given by a simple condition of R IB < 1 based on the buoyancy Reynolds number R IB .
The effect of this transition is also evident in Lagrangian statistics. For acceleration statistics, the effect of imposed RaS appears very weak before the transition, suggesting the imposed anisotropy does not affect the small-scale isotropy significantly. The probability density functions (PDFs) of acceleration in both horizontal and vertical direction display wide tails very similar to isotropic turbulence, and the variance of acceleration can be approximately scaled with Kolmogorov variables. The Lagrangian autocorrelations and frequency spectra also appear to behave in a similar fashion. However, the sharp transition at F r 0.07 leads to a very strong anisotropy between the horizontal and vertical directions, with the particle motion strongly dominated by RaS instead of turbulent eddies. The PDFs of acceleration show suppressed tails, and acceleration variance shows no clear scaling. Moreover, the autocorrelations and frequency spectra clearly reflect the dominant roles of frequencies N and f .
On the other hand, Lagrangian velocity statistics are always affected by the imposed RaS. Anisotropy is already evident for weak RaS; however the degree of anisotropy becomes very strong with the said transition at F r 0.07. This is most readily seen in integral time scale based on the autocorrelations from vertical velocity and hence also reflected in long time behavior of meansquare displacement of particles in the vertical direction. Similar effects had been observed in the case of purely stratified flows [19]. In such flows, it is known that stratification leads to formation of horizontal layers, that inhibits the transport in the vertical direction. In the presence of both RaS, slanted layers appear which inhibit, to some degree, transport in both the vertical and the horizontal directions.
It is interesting to consider that the transitional behavior observed in this work is possibly also linked to the enhanced vertical velocities that develop in strongly stably stratified turbulent flows, such as in the nocturnal planetary boundary layer [53], and warrant further examination. The resulting dispersion of Lagrangian particles must feel these hot spots of strong vertical velocity and strong local dissipation linked to localized shear instabilities, as also observed in the ocean. Several studies point to a local Richardson number, measuring stratification with respect to the internal shear associated with vertically sheared horizontal winds, to be mostly in the vicinity of the value for linear instability [28,[54][55][56]. This leads to flows which are altogether close to criticality, and strongly anisotropic, intermittent and dissipative [55]. We expect that the observed similarity between the statistical properties of acceleration when the flow is dominated by eddies (for R IB > 1), and those in HIT, may be very useful in modeling the dispersion of particles in turbulent flows in the presence of RaS.
We conclude by mentioning that many interesting issues remain to be addressed for particle motion in such flows. For example, the problem of separation of pairs, or of clusters of particles, both forward and backward in time [57,58] or considering the effects of particle inertia or molecular diffusion. In addition, exploring more values of N/f and P r at possibly higher Reynolds numbers will help in better understanding the parameter space and how the results from current work translate to practical situations [25,45].

ACKNOWLEDGMENTS
The computing resources utilized in this work were provided by PSMN at the Ecole Normale Superieure de Lyon. D. Buaria and A. Pumir acknowledge support from the European High-