Shape and size of large-scale vortices : a universal fluid pattern in geophysical fluid dynamics

Planetary rotation organizes fluid motions into coherent, long-lived swirls, known as large scale vortices (LSVs). LSVs are ubiquitous in nature, and their shape and size are expected to control their effect on the dynamics and long-term evolution of geophysical and astrophysical fluids. By using high-resolution direct numerical simulations, here we show for the first time that the shape of LSVs in rapidly-rotating mixed convective and stably-stratified fluids, which approximates the two-layer, turbulent-stratified dynamics of many geophysical and astrophysical fluids, is universal and that their size can be predicted. Specifically, we show that LSVs emerge in the convection zone from upscale energy transfers and decay as they penetrate into the stratified layer by thermal wind balance, thus taking the shape of a depth-invariant cylinder in the turbulent layer and of a penetrating half dome in the stable one. Furthermore, we demonstrate that when LSVs penetrate all the way through a stratified layer bounded by a solid boundary, they saturate by boundary friction. We provide a prediction for the penetration depth and maximum radius of LSVs as a function of the LSV vorticity, the stable layer depth and the stratification. Our results suggest that while turbulent vortices can penetrate far into the stratified layers of atmospheres and oceans, they should stay confined within the convective layers of Earth's liquid core and of the Sun.


I. INTRODUCTION
Large-scale vortices (LSVs) are a key component of geophysical and astrophysical fluids. They are generated by a myriad of processes, ranging from the instability of currents and fronts in oceans [1] to tropical cyclogenesis in the atmosphere [2]. In oceans, LSVs have O(1 − 100)km diameter, weeks to years lifespan [3,4], and they can transport ocean mass, heat and CO 2 over long horizontal [5][6][7] and vertical [8] distances. They also influence the background flow [9] and significantly affect plankton productivity and chlorophyll distribution in surface waters [10,11]. Planetary atmospheres showcase a wide range of LSVs [12], including long-lived large planetary-scale vortices that control the global circulation and climate (e.g. polar vortex on Earth and Jupiter's Great Red Spot) as well as smaller cyclones with O(100)km diameter on Earth [13] that can have devastating consequences. Earth's outer core, which is made of turbulent liquid iron that powers the Earth dynamo, is also expected to feature numerous LSVs of various sizes [14], such as the high-latitude geomagnetic flux patches [15], as well as a large-scale north polar vortex [16]. LSVs are also found in the solar photosphere [17], and are expected to exist in accretion disks [18] and potentially play an important role in planet formation [19].
LSVs typically result from the breakup of large-scale flows or from upscale energy transfers that feed on small-scale waves and turbulence. In shallow fluid layers that are considered two-dimensional, an inverse cascade guarantees a flux of energy from small scales to LSVs [20,21]. However, stars and Earth's outer core can hardly be considered shallow, and LSVs in Earth's atmosphere and oceans often have complicated vertical structures, such that three-dimensional theories are required for realistic predictions. In recent years, it has been shown that rapid rotation enables upscale energy transfers in fully three-dimensional turbulent convection, with a barotropic large vortical mode emerging from the turbulent eddy field [15,[22][23][24]. The convection must be turbulent but also strongly constrained by rotation for the LSV to emerge, a regime known as geostrophic turbulence. If geostrophic turbulence is common in geophysical and astrophysical fluids, dedicated simulations could unravel the fundamental characteristics of many of the LSVs in nature. Specifically, the relationship between core pressure anomaly, maximum velocity and size could be investigated rigorously, and help predict the impact of LSV not only in the oceanic and atmospheric contexts [25], but also in planetary physics and astrophysics. We remark that previous works have focused on simulations of fully-convective fluids and used free-slip boundaries: in this context, no physical process has been found to saturate the growth of LSVs (i.e. LSVs always reach the box size), except for magnetism in a recent study [26].
We extend previous studies of LSVs in fully-convective fluid systems to LSVs in fluids that are self-organized in a turbulent layer next to a stably-stratified fluid region. Our aim is to investigate the shape and size of a generic model of LSVs similar to eddies in the surface ocean mixed layer penetrating into the thermocline, of cyclones in the Earth's turbulent planetary boundary layer reaching into the upper troposphere and stratosphere, and of LSVs in the convection zone of stars and planetary liquid cores overshooting in adjacent stable (radiative) layers. In the Earth's core context, evidence of a stably-stratified layer at the core-mantle boundary [27,28], or adjacent to the inner core [29,30], is recent and has prompted significant interests owing to its potential influence on the geodynamo [31] and core flows [32]. Past studies of penetrating vortices in mixed convective-stably-stratified cores, e.g. [33][34][35], are limited to a regime dominated by bulk viscosity that is unlikely to be relevant to planetary dynamics [15]. Our results may also be applicable for subsurface oceans, e.g. on Enceladus, and subglacial lakes, when a stratified layer exists at the bed-water or ice-water boundary, as may be the case for subglacial lakes in Antarctica close to the surface [36].
Here we show how finite stable fluid layers and boundary friction can control the maximum size of LSVs. This is a result of significant importance since the saturation of upscale energy transfers is not universal but depends on the dissipative or dispersive mechanisms at play at large scales. We demonstrate that the key features of LSVs, including core pressure anomaly, LSV diameter and maximum azimuthal velocity, satisfy thermal wind balance. We show that the stability and depth of the stable layer control the diameter and extent of penetration of LSVs into the stable layer, which ultimately sets the LSVs' potential to promote vertical exchanges across density interfaces. These two effects are investigated systematically using a suite of direct numerical simulations (DNS) of the Navier Stokes equations with high resolution and long integration time. We derive an aspect ratio for the penetrating, stably-stratified part of LSVs and deduce an approximate penetration depth and maximum size of LSVs in different geophysical and astrophysical contexts.

II. PROBLEM FORMULATION
We investigate the dynamics of a local fluid domain confined between solid top and bottom boundaries and rotating at constant (Coriolis) frequency f about the vertical z axis (ẑ the upward pointing unit vector). The governing equations are the Navier-Stokes equations under the Boussinesq approximation and we assume horizontal periodicity. Bottom boundary conditions are always imposed temperature (T = 1 in dimensionless space) and free-slip velocity, while top boundary conditions are imposed temperature (T = T t < 0) and either free-slip or no-slip velocity. We use a nonlinear equation of state (cf. equation (1d) below) in order to obtain a self-organized, mixed convective and stably-stratified fluid. Similar to water near its density maximum at 4 • C, the thermal expansion coefficient is temperature-dependent and changes sign [37]. Specifically, it is piecewise-constant, positive for T > 0 and negative for T < 0. Thus, the density is maximum at T = 0 away from the bottom boundary and the lower layer of fluid is convectively unstable whereas the top layer is stably stratified [38,39]. We use the height of the convective layer and the thermal diffusive time as reference length and time scales. The dimensionless equations for velocity u = (u, v, w), pressure p, temperature T and density anomaly ρ are in a Cartesian (x, y, z) frame of reference where H is the Heaviside function, P r = ν/κ is the Prandtl number, Ra = α s g∆d 3 /(νκ) is the Rayleigh number, Ek = ν/(f d 2 ) is the Ekman number and S is the stiffness parameter, with ν the viscosity, α s the thermal expansion coefficient for the convecting fluid, g the gravity and ∆ the temperature difference driving the convection. We note b = −P rRaρ the buoyancy. Dimensional variables can be recovered from the dimensionless ones using d, d 2 /κ, ∆ and ρ 0 α s ∆ as characteristic length, time, temperature and density scales, with ρ 0 the reference density of the fluid. The control parameters are the horizontal size of the box L = L x = L y (in units of convective layer depth), Ra, Ek, P r, the dimensionless stable layer depth H and the background buoyancy frequency N = −ST t /H (also known as Brunt-Väisälä frequency) with S the ratio of the thermal expansion coefficient in the stable layer to the thermal expansion in the convective layer (also known as stiffness parameter). Here, we set L = 4 and we select P r = 1, Ra = 2 × 10 8 and Ek = 10 −5 such that the lower convective layer is in the regime of geostrophic turbulence of fully-convective fluids which features LSVs [40]. For each L z considered, we adjust T t such that the convection zone extends from z = 0 to z ≈ 1 [39]. Thus, we write L z = 1 + H with H ≥ 0 the stable layer thickness. We use the high-performance, opensource pseudo-spectral simulation code DEDALUS [41] in order to solve the governing equations. We use Chebyshev and Fourier modes in the z and x, y directions, respectively, and a 2-step implicit/explicit Runge-Kutta scheme for time integration. The CFL condition is set to 0.45 and the time step is typically O(10 −7 ) (cf. resolution in table I).
We investigate the effect of the buoyancy frequency N (proxy for stratification strength) on the penetration of turbulent LSVs into the stable layer and we demonstrate that LSVs saturate when they penetrate through the entire stable layer depth and reach the top no-slip boundary. We run the simulations long enough such that all results presented are at quasi steady-state, i.e. the heat flux is constant throughout the depth of the whole fluid and LSV properties are at statistical equilibrium. Key parameters of the simulations are presented in table I, and additional figures are given in the supplementary information (SI). Note that all LSVs that are coherent in our simulations are cyclonic, in agreement with previous DNS of fully-convective rotating fluids under Boussinesq approximation [23].

A. Importance of stably-stratified layers
We first show in figure 1 three-dimensional snapshots of the horizontal velocity amplitude V = √ u 2 + v 2 at steady state. In fully-convective simulations without a stable layer, a LSV emerges when the top boundary is free-slip (figure 1A * ), but not when the top boundary is no-slip (cf. figure 1A). Thus, boundary friction inhibits upscale energy transfers in fully-convective fluids, to the point that, as shown by previous studies [42], large-scale barotropic vortices cannot be obtained in current DNS (i.e. which are limited to relatively high viscosity) with no-slip boundaries. With a stable layer (H = 0), we find that one or several LSVs always emerge for the same convective parameters as in figure 1A, even with a no-slip top boundary (cf. one LSV in figure 1B and several smaller LSVs in figure 1C). This means that stable layers protect upscale energy transfers and LSVs against boundary friction, which is a fundamental and important result for planetary cores and potentially for Earth's oceans and subsurfaces oceans. It shows that subadiabatic layers of planetary cores and oceans' pycnoclines can play an important role in the dynamics of LSVs by protecting them against boundary friction at e.g. the core-mantle boundary or the seabed. It may be noted that LSVs are expected to be robust against no-slip boundaries in reduced models of fully-convective fluids assuming asymptotically-large rotation and turbulence intensity [43]. Therefore, a stable layer tampering boundary friction may not be always necessary, but it still broadens the domain of existence of LSVs to cases accessible to DNS and possibly laboratory experiments [44]. We recall that the bottom boundary is free-slip in all our simulations, since LSVs cannot emerge in a convective fluid directly adjacent to a no-slip bottom boundary for our choice of parameters.

B. Horizontal saturation
LSVs in nature grow to a finite size, i.e. saturate, either because there is a physical mechanism that prevents their growth beyond a certain point or because they reach the boundaries of the geophysical or astrophysical fluid domain. Previous studies of fully-convective Cartesian fluid domains with free-slip boundaries have always reported LSVs growing to the box size [23,24]. This is a severe limitation to the application of existing numerical local models to natural cases, since the box size in periodic simulations is not a real physical quantity. Here, we demonstrate that boundary friction through a stably-stratified layer provides a natural saturation mechanism for LSVs, and that the final natural diameter depends on the stratification strength N and depth H of the stable layer. Figures 1B,C clearly show the sensitivity of the natural diameter of LSVs with N (all other parameters being the same). In figure 1B the stratification is strong and the LSV fills up the entire domain, suggesting that the natural LSV diameter is large, larger in fact than the horizontal extent of the domain. In figure 1C, on the other hand, the stratification is weak and several LSVs co-exist, merge and split but on average do not grow bigger than about a third of the domain size, suggesting that the LSVs saturate naturally at a moderate diameter and do not experience numerical confinement.
In order to assess which simulations feature domain-filling LSVs (i.e. confined numerically) and which simulations feature LSVs saturating naturally, we show in figure 2A the integral length scale L 0 , a proxy for LSV diameter, in the middle of the convection zone. The integral length scale is given by where k is the horizontal wavenumber and a hat denotes Fourier transform in (x, y). In all cases, we see that L 0 has saturated by t = 0.1 (recall that t is normalized by the long thermal diffusive time scale). The fully-convective simulation C has L 0 ≈ 0.5 (solid black line), which is much smaller than L 0 in all other cases. This is because there is no LSV in C due to the no-slip condition, as seen in figure 1A. Conversely, for simulation C * (dashed black line), L 0 ≈ 3.3 and is roughly the maximum attainable since in this case the cyclone saturates at the size of the numerical domain L = 4. With a stable layer and a no-slip top boundary (colored solid lines), we find that L 0 at steady-state increases with the stratification N (blue to orange to green) and with the stable layer depth H (thin to thick lines).
The simulations with the strongest stratification (S 0.5 and S 1 ) and with the thickest stable layer (M 2 ) have L 0 ≈ 3.3, i.e. feature a unique LSV that has reached the domain size. The effect of the stable layer depth on the number and diameter of LSVs in simulations with moderate stratification can be seen in figure 2B-D where we show the horizontal velocity V in the middle of the convection zone at t ≈ 0.1: clearly, LSVs saturate at smaller sizes when the stable layer becomes shallower (figure 2B to 2D). It is worth noting that LSVs saturate naturally only when the top boundary is no-slip and provides friction. Indeed, while L 0 ≈ 1.2 for W 0.5 with no-slip, L 0 ≈ 3.3 for W * 0.5 with free-slip and the LSV fills up the entire domain (i.e. saturate numerically). In fully-convective systems it has been shown that the box-filling LSVs can be replaced with large-scale jets when the domain aspect ratio is changed [45,46]. Our results suggest that moderate-size, penetrating LSVs should be robust against such changes since they saturate at diameters smaller than the horizontal extent of the numerical domain.

C. Shape of penetrating LSVs
In order to understand why weak stratification and small stable layer depth (resp. strong stratification and large layer depth) lead to small (large) LSV diameters, we now provide a phenomenological description of the axisymmetric shape of LSVs. Importantly, we remark that the thermal wind balance is satisfied, at least at first order, by all LSVs in our simulations, such that, in a cylindrical coordinates system centred on a LSV, we have, after time and azimuthal averaging, where v θ is the azimuthal velocity and b is the buoyancy anomaly, i.e. b = b − b ∞ with b ∞ the buoyancy in the far field (cf. Appendix A). Note that the vertical vorticity ζ is related to the velocity via rζ = ∂ r (v θ r), such that ∂ z ζ and ∂ z v θ have generally the same sign in a LSV. We show in figure 3A a schematic of a LSV in mixed convective and stably-stratified fluid, which is based on the time-averaged and azimuthally-averaged map of normalized vertical vorticity obtained in DNS (cf. figures 3B,C). We note the radius of maximum azimuthal velocity at the base of the stratified vortex cap, h the e-folding decay height of azimuthal velocity in the stratified layer (penetration depth for short) and δ the restratification depth. The black dashed line shows the interface between the convection zone and the stably-stratified layer. We represent the LSV in figure 3A as a cylinder of depth-invariant vorticity within the convective layer and as a half dome of vertically-decaying vorticity in the stably-stratified fluid, which we call the stratified vortex cap. The large vorticity inside the LSV inhibits turbulence compared to the outside in the convective layer [47], such that there is less and less mixing toward the LSV centre. This results in a vertical temperature gradient steepest at the LSV centre [24], and, accordingly, a downward depression of the isothermal of maximum density (black dashed line) also toward the LSV centre. We call the extent δ by which the stratified vortex cap sinks into the convective zone the restratification depth, in reference to the restratification of the oceanic surface layer due to eddies [48]. The decrease (in magnitude) of the vertical temperature gradient with radius results in a negative temperature anomaly, T = T − T ∞ , in the LSV centre. This anomaly is shown by the light red-coloured cone in figure 3A and is small, as is the buoyancy anomaly b = T < 0, in most of the convective layer. As a result, the LSV roughly satisfies the Taylor-Proudman theorem, i.e. is depth-invariant, in the convective layer (cf. equation (3)). The negative temperature anomaly increases with height, such that at and above the base of the stably-stratified layer, it translates into a positive and potentially large buoyancy anomaly b = −ST . This positive buoyancy anomaly drives the decay of the azimuthal velocity with height above the black dashed line according to the thermal wind balance (equation (3)), which is why the stratified LSV has a half-dome shape. When S increases, i.e. the stratification becomes stronger, b increases, such that the aspect ratio h/ of a LSV must decrease in order to satisfy the thermal wind balance. This explains why in a strongly-stratified fluid LSVs appear as wide weakly-penetrating columns (cf. figure 1B), while in a weakly-stratified fluid they appear as tall narrow columns (figure 1C).
Figures 3B,C show the vertical vorticity for simulations M 2 (figure 3B) and M 0.5 (figure 3C), i.e. which have a deep and shallow stratified layer, respectively, but same parameters otherwise. As described above, the stratified vortex cap has a positive buoyancy anomaly in both cases (as shown by the gray contours), which is balanced by a vorticity decay with height above the convective-stable interface (dashed line). However, while the penetration of the vortex cap is small compared to the stable layer thickness H in figure 3B, the penetration is large enough compared to H in figure 3C such that the LSV is confined vertically. The maximum vorticity does not change significantly between the two simulations and the buoyancy anomaly is smaller in figure 3C than in figure 3B (cf. in-line numbers). Thus, |∂ z v θ | is larger for a vertically-confined LSV than for a vertically-unconfined LSV, which means that confined LSVs must decrease in diameter (compared to their unconfined counterparts), i.e. such that |∂ r b | increases, in order to maintain thermal wind balance. As a result, boundary friction makes the LSVs saturate naturally in general and in particular in figure 3C, because it imposes a sharp vorticity decay that can only be balanced by a reduction of the LSV diameter. It can be noted that the horizontal narrowing of vertically-confined LSVs does not apply when the top boundary is free-slip since in this case the vorticity doesn't decay any quicker than when it is unconfined.  The aspect ratio of the stratified vortex cap, α = h/ , is a function of the normalized stratification strength N /f , with N = f EkN/P r the dimensional buoyancy frequency, and the Rossby number of the LSV, i.e. Ro = Ek(v 0 θ /P r)/ with v 0 θ the maximum azimuthal velocity at the base of the stratified vortex cap. An approximate expression for α(Ro, N /f ) can be derived from the hydrostatic and cyclo-geostrophic equations, which are slightly more relevant in our case of small but finite Rossby than the thermal wind balance (3) and which read with p the pressure anomaly (again using a cylindrical coordinates system centred on a LSV core and time-and azimuthally-averaged variables). The relationship between α, Ro and N /f arises from the requirement that the pressure anomaly in the core due to the cyclo-geostrophic flow (cf. equation (4a)) must be the same as the pressure anomaly due to the positive buoyancy anomaly of the stratified vortex cap (cf. equation (4b)). This is a type of consistency condition that leads to an expression for α(Ro, N /f ) that depends on the radial profile of azimuthal velocity and on the vertical profile of buoyancy. The formula for α(Ro, N /f ) was previously derived for vortices in fully-stratified fluids [49] and lenticular vortices at the ocean surface [50], and here we derive it for turbulent LSVs penetrating in a stably-stratified fluid. We find that the radial profile of LSVs in DNS matches reasonably well with the radial profile of shielded monopoles [51], and that the vertical profile of buoyancy anomaly is well approximated by a constant substratified bottom with an exponentially-decaying cap. This yields the formula (cf. appendix B) with a 1 and a 2 parameters of order unity, given by with Γ the Gamma function, µ the steepness parameter of the velocity profile in r, b 0 the buoyancy anomaly at the base of the stratified vortex cap and N 0 < N the stratification strength of the stratified vortex cap inside the convective layer. Equation (5) is derived under the assumption that the LSV is in an infinitely deep and wide stably-stratified fluid, i.e. such that the stratified vortex cap doesn't reach the top boundary. In our simulations, we have both vertically confined and unconfined LSVs and our numerical box has a finite horizontal extent, such that equation (5) cannot be expected to be satisfied exactly. Nevertheless, we show in figure 4A that the aspect ratio h/ measured directly from the velocity profile v θ in DNS matches very well with the theoretical prediction (5) based on the problem parameters, such that the formula is applicable for both unconfined and confined LSVs.
From equation (5) we can obtain an approximate expression for the maximum diameter of LSVs penetrating in a stably-stratified fluid. The maximum diameter of LSVs is the diameter of LSVs that are confined vertically and saturate naturally by boundary friction (in the absence of other saturating mechanisms). We show in figure 4B the radius of LSVs in our simulations as a function of H/α with α given by equation (5). To the left of the diagram, i.e. where H is small, we have the results of LSVs that are confined vertically and saturate naturally. For such LSVs, we find that ≈ H/(2α) = max , which we therefore define as the maximum radius of LSVs. To the right of the diagram, where H is large and LSVs saturate horizontally at the box size before they reach the top boundary, we find that < max , as expected. Note that based on figure 4B is in fact close to max + with ≈ 0.18 a small correction in the limit H/α → 0, which may be due to complicated boundary layer effects that are neglected in the present work.
We find that µ ≈ 1 such that a 2 ≈ 2/3 in all simulations, i.e. for both vertically confined and unconfined LSVs, and that a 1 ≈ 5/2 for unconfined LSVs but varies with the problem parameters, i.e. a 1 ∈ [2, 5], for confined LSVs (cf. appendix B). Thus, in the Discussion section we will use for the penetration depth of unconfined LSVs and for the maximum radius of confined LSVs the approximate formulas i.e. such that max is an upper bound based on our DNS results (i.e. using the minimum of a 1 ).

IV. DISCUSSION
We have shown that upscale energy transfers in geostrophic turbulence generate LSVs in mixed convective and stably-stratified fluids. The LSVs are depth-invariant in the convective layer and decay in the stable layer by thermal wind balance because the LSV core is positively buoyant. The growth of LSVs stops when LSVs penetrate through the entire stable layer depth and reach the top no-slip boundary. Thus, in addition to the well-known beta-effect [e.g. 52] and to the presence of strong magnetic field [26], boundary friction across a stably stratified layer constitutes a physically relevant saturation mechanism to quench the inverse cascade of rapidly rotating, convective turbulence in natural systems.
LSVs studied in this work may be considered a simplified model of cyclones in Earth's atmosphere [53], and in particular of warm-core tropospheric cyclones penetrating into the stratosphere [54], of eddies in Earth's oceans [55,56], and of LSVs in Earth's outer core [14] and stars. Earth's atmosphere, oceans, outer core and stars have different fluid properties, such that the aspect ratio of the stratified cap of LSVs, and the penetration depth, depend on the geophysical or astrophysical fluid of interest. We give in figure 5 different values of the aspect ratio α 0 of equations (7) in (N /f, Ro) plane, and we highlight regions relevant to Earth's atmosphere, oceans, outer core and stars. The atmosphere and oceans are relatively strongly stratified (i.e. N /f ≥ 1) and have high Ro and moderate Ro, respectively. On the other hand, Earth's outer core is expected to be moderately stratified with small Ro, and stars have weak stratification and moderate Ro. As a result, LSVs are expected to be wide and weakly-penetrating in Earth's outer core and stars, while moderately-penetrating in Earth's atmosphere and oceans.
For LSVs in Earth's atmosphere, if we take Ro ∼ 1 and ∼ 100km, equation (7) yields h ∈ [2.5, 25]km for N /f ∈ [10,100]. Thus, our model predicts that atmospheric LSVs can reach far into the stratosphere, and potentially all the way to the ozone layer found at ≈ 25km when the turbulent planetary boundary layer is deep and atmospheric stability is low. In Earth's oceans, mesoscale eddies have typically ∼ 100km and Ro ∼ 10 −2 (based on rms velocity ∼ 10cm/s) [57], and submesoscale eddies have typically ∼ 10km and Ro ∼ 1 [58]. Thus, the penetration depth of both eddy types is the same, i.e. h ∈ [0. 25,25]km for N /f ∈ [1, 100], which shows that surface eddies can penetrate relatively far into the thermocline and potentially reach the seabed, especially in weakly-stratified waters on the continental shelf. In the Earth's core, if we take ∼ 30km and Ro ∼ 10 −6 for the diameter and Rossby number of the most intense LSVs, as suggested in a recent study [14], we find h ∼ 3m for N /f = 1, which is a typical value used in previous works, e.g. [27]. This result suggests that upwellings and downwellings inside dominant LSVs in Earth's core do not promote exchanges of chemical species between the convection zone and far into the stably-stratified layer, unlike LSVs in the atmosphere and oceans. We note that non-dominant LSVs in Earth's core may penetrate farther into the stably-stratified layer. Previous works on fully-turbulent outer core dynamics studied LSVs at both planetary scale ∼ 1000km with Ro ∼ 10 −5 and smaller scales ∼ 100km with Ro ∼ 10 −4 [59,60]. For such LSVs and N /f = 1, we find h ∼ 10km and h ∼ 2.5km, respectively. In stars, Ro ∼ 1 is relevant for supergranulation [61] and N /f ∼ 10 3 is a reasonable estimate for the stratification of the Sun [62]. If ∼ 0.1R * with R * the star radius, then h ∼ 10 −4 R * . Thus, LSVs in stars similar to the Sun are weakly-penetrating and cannot go through the tachocline, which is on the order of one percent of the stellar radius for the Sun [63].
In Earth's oceans the thermocline protect LSVs from the seabed, and in Earth's outer core stably-stratified layers may protect LSVs at both the inner-core and core-mantle boundary. The seabed and solid boundaries around Earth's outer core provide friction, which may play a role in the saturation of LSVs. The thermocline of Earth's oceans is on the order of a few km deep, H ∈ [1, 10]km, which means that the maximum diameter of LSVs in Earth's oceans, for a moderate stratification of N /f = 10, is max ∈ [25,250]km for mesoscale eddies (Ro = 10 −2 ) and max ∈ [2.5, 25] for submesoscale eddies (Ro = 1; cf. equation (7)). Since the lower bound of max lies in the range of observed ocean eddies, our work predicts that the seabed may play a role in limiting the size of ocean LSVs. The thickness of the stably-stratified layers in Earth's core is poorly constrained. Recent studies use H ∼ 100km or more [30,31]. For H ∼ 100km, we find, for N /f ∼ 1, max = 25000km for Ro = 10 −6 and max = 2500km for Ro = 10 −4 . max is larger than the radial extent of Earth's outer core for both low and high Ro, such that boundary friction is unlikely to be the relevant saturation mechanism for LSVs inside the Earth. This implies that for studies discarding the stable layer, a stress-free boundary condition for the turbulent flows may be more appropriate than a no-slip condition.
Finally, we note that our model neglects compressibility effects (which may be important for LSVs in the atmosphere, outer core and stars; cf. [64,65]), radiative transfers (atmosphere and stars), moist dynamics (atmosphere) and magnetic field (outer core and stars). Our work also neglects the dynamics of the Ekman layer at the top of the stably-stratified fluid, which may yet affect the prediction for the shape and size of vertically-confined LSVs, as suggested by the variability of a 1 in equation (5). Future investigations taking into consideration one or several of the effects neglected in this study will help further our understanding of LSVs in nature.
with v 0 θ (z) the maximum azimuthal velocity, (z) the radius where the velocity is maximum, and µ(z) ∼ O(1) the profile steepness. We measure v 0 θ and from the DNS results at each z and obtain µ(z) by least-square fit for r ∈ [0, r max ]. We show µ in figure 3 in SI: µ is roughly equal to unity in the convection zone (exponential decay of vorticity in r) and then increases to approximately two in the stable fluid (Gaussian decay, which is typical of vortices in stably-stratified fluids, e.g. [49]). We find that there is some variability of µ depending on r max (i.e. the extent over which we perform the best fit), which is not surprising since in our simulations the LSVs cannot relax to infinity but are instead horizontally periodic. We take the average of the best-fit values for µ for r max ∈ [1, 1.5].
The pressure anomaly p satisfies the geostrophic equation (A1a) and p = 0 at r = z = ∞, i.e. assuming an unbounded fluid domain. The far-field pressure is p ∞ = p − p and we obtain the far-field buoyancy b ∞ from equation (A1b) with p substituted by p ∞ . Note that we require p ∞ and b ∞ to be piecewise second-and first-order polynomials (the buoyancy must have a purely diffusive profile) joined at z = z ∞ with z ∞ ≈ 1 the turbulent-stable interface position in the far field. In the results presented z ∞ is let arbitrary and obtained by best fit but it can be taken equal to one without substantial changes.
In our simulations we find that the buoyancy anomaly is well approximated along the rotation axis by a profile of the form with b 0 the buoyancy anomaly at z = z 0 , N 0 the density restratification due to the LSV in [z 0 , z 1 ], z 1 ≈ 1 the profile transition height, and h the overshooting depth parameter; z 0 and b 0 are taken from the simulation results, while N 0 , z 1 and h are obtained by best fit for z ≥ z 0 . Note that setting z 1 = 1 adversely impacts the fit for the cases with weak stratification. The depth of restratification of the fluid below z 1 is noted δ in figure 3A of the main text. The velocity v θ (r, z 0 ) and its fit (A2) as well as the buoyancy anomaly b (0, z) and its fit (A3) match well as shown in figures 4 and 5 of the SI. with Γ(·) the Gamma function and Ro = Ek(v 0 θ /P r ) the Rossby number based on v 0 θ and at the base of the stable layer (P r appears because we use the thermal diffusive time scale for normalization). Integrating (A1b) for the pressure and buoyancy fit (A3) vertically along r = 0 then yields Assuming that the pressure anomaly far from the vortex is 0, i.e. p (∞, z δ ) = p (0, ∞) = 0, and equating (B1) and (B2) we finally obtain for the aspect ratio squared which yields (5) of the main text with EkN/P r rewritten as N /f (N the dimensional buoyancy frequency), and the parameters a 1 and a 2 are given by We find that a 1 ≈ 5/2 is approximately constant for unconfined LSVs, while for confined LSVs a 1 ∈ [2, 5] shows some variability (figure 6 in SI). Furthermore we find that a 2 ≈ 2/3 in all simulations (also figure 6 in SI). The variability of a 1 for confined LSVs comes from the fact that (i) the fit proposed for the buoyancy anomaly does not accomodate