Mechanisms for Limiting the Depth of Zonal Winds in the Gas Giant Planets

We solve the above equations in a spherical shell with inner radius r_i = 4 and outer radius r_o = 5 in the Boussinesq limit, i.e., $\tilde{\rho }$ = 1. Spherical coordinates (r, θ, ϕ) are used and only axisymmetric flow is considered, i.e., ∂/∂ϕ = 0. The dimensionless background codensity gradient ${C}_{o}^{{\prime} }$ equals one at r_i and decreases mildly with r so that the total flux ${r}^{2}{C}_{o}^{{\prime} }$ remains constant up to r − r_i < 0.32. For r − r_i > 0.48, it is zero with a smooth transition in between. The diffusivity κ is set to the reference value κ_o in the upper part of the shell (i.e., nondimensional value is one), and is reduced tenfold in the stable lower region, also with a smooth transition. The driving force is specified as ${{\boldsymbol{F}}}_{{\boldsymbol{D}}}={A}_{f}\,f(r)\gamma (s)h(s){{\boldsymbol{e}}}_{\phi }$, where A_f is the amplitude of forcing, ${{\boldsymbol{e}}}_{\phi }$ is the unit vector in azimuthal direction, s = r sin(θ) is the distance to the rotation axis (or cylindrical radius). The forcing is restricted to radii larger than r_f, which is set here to 4.7, i.e., forcing is limited to the top 30% of the shell. Setting $\bar{r}=(r-{r}_{f})/(1-{r}_{f})$, the radial dependence is, for r > r_f:

$$\begin{eqnarray}&&f(\bar{r})=3{\bar{r}}^{2}-2{\bar{r}}^{3}.\end{eqnarray} \tag{ 10 }$$

Figure 1(a) shows the radial variation of the background codensity, of the codensity diffusivity, and of the forcing function. Seven bands of alternating zonal flow inside the tangent cylinder $s\lt {r}_{i}$ are driven in each hemisphere by specifying an s-dependence of the forcing as

$$\begin{eqnarray}&&\gamma (s)=\sin (7.5\pi \,s/{r}_{o}).\end{eqnarray} \tag{ 11 }$$

Friction outside the tangent cylinder is less than it is inside. To avoid excessive differences in zonal flow velocities between the two regions, the forcing function is attenuated in the s-direction, setting h(s) to one for s < 3, to $1-0.5{\left(s-3\right)}^{2}$ for 3 < s < 4, and to $0.5{\left(5-s\right)}^{2}$ for s > 4. The complete pattern of the forcing function is displayed in Figure 2(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The boundary conditions are impermeable and shear-stress free at r_o and at r_i. We note that the forcing function implies a net torque on the shell; in the absence of magnetic coupling to the interior, which is immobile in the rotating reference frame, this would imply a continuous spin-up. That is taken care of by simply removing the solid-body rotation part from the solution. The condition for the codensity perturbation is c = 0 at r_o and dc/dr = 0 ar r_i. The magnetic field at r_o is matched with a potential field above the shell. An axial dipole field with amplitude B_o at the poles is imposed as a boundary condition at r_i. All other magnetic field components are matched with a field that is obtained by solving Equation (2) with ${\boldsymbol{u}}=0$ and η = 1 in the region inside the sphere. This means that the electrical currents associated with a toroidal magnetic field, which is induced by the flow in the shell, can close through the highly conducting interior. The discontinuity in horizontal velocity at r_i between the flow in the shell and the immobile core implies a complex nonlinear jump condition for the radial derivative of the toroidal field at the interface; see Equations (5.8)–(5.10) in Kuang & Bloxham (1999).

To solve the equations numerically, a pseudo-spectral dynamo code that expands all variables in spherical harmonic functions in the angular variables and in Chebychev polynomials in radial direction is used (Christensen & Wicht 2015). Here, it is stripped down to only axisymmetric modes. We use 341 Legendre polynomials in the θ-direction and 97 Chebychev polynomials in the radial direction. By using a grid-stretching technique, we condense the radial grid points in the transition region between the stable and the unstable parts. The system is marched in time until a stationary solution is reached.

In all models presented here, the Ekman number has been set to E_κ = 10⁻⁵, and the Prandtl number Pr and the Roberts number (at the bottom) q are one. The magnetic diffusivity increases exponentially

$$\begin{eqnarray}&&\eta =\exp (\mathrm{ln}({\rm{\Delta }}\eta )(r-{r}_{i})),\end{eqnarray} \tag{ 12 }$$

where the change from the inner boundary to the outer boundary is by a factor ${\rm{\Delta }}\eta ={10}^{8}$. The setup and the choice of parameters was not made with the aim to represent conditions for Jupiter or Saturn as closely as possible. Instead, the models are meant to explore and illustrate general principles that govern the depth-dependence of zonal flow.

Figure 2(a) shows the velocity field in the nonmagnetic case without a stable layer. The zonal flow satisfies the Proudman–Taylor theorem almost perfectly, i.e., it hardly varies with z. The meridional circulation is weaker than the zonal flow. In the top region, the force driving the zonal flow is balanced by the Coriolis force associated with the meridional circulation according to Equation (13), hence it implies a flow in the s-direction. At depth, the direction of u_s must reverse because of mass conservation, and here the Coriolis force must be balanced by viscous friction in the bulk of the fluid, resulting in a distributed return flow, when Lorentz forces are unavailable.

We next consider nonmagnetic case 2 with a strongly stable layer ($\tilde{N}\approx 4$) in the lower part of the shell. Inside the tangent cylinder, the jet flow in the neutrally stable region is slightly stronger than in the absence of a stable region (Figure 2(c)). In the stable layer, it decreases with depth in amplitude, although it remains significant down to the lower boundary. The meridional circulation grazes the top of the stable region where it creates a codensity anomaly by lateral advection, which is strongest in the upper part of the stable layer (Figure 2(d)). The zonal flow velocity decreases with depth, due to the associated thermal wind effect given by the first term on the rhs of Equation (14). The jets become slanted, i.e., they extend in the radial direction in the stable layer rather than along the z-axis. The ratio of the rms velocity on the bottom boundary to that at the surface inside the tangent cylinder (s < 3.95) is listed in Table 1 under the entry R_u. In the cases with a stable layer but without magnetic field, R_u is still approximately 0.5 at $\tilde{N}\approx 4$. This simulation suggests that, in a real planet, a stable layer alone—without some additional resistive effect on the flow—may be rather inefficient in reducing zonal flow velocities with depth.

We now turn to cases 3 and 4, which lack a stable layer but have an exponentially varying electrical conductivity and an imposed dipolar field. While the jet outside the tangent cylinder remains strong, the amplitude of the zonal flow inside the tangent cylinder is reduced by more than two orders of magnitude compared to cases 1 and 2 (Table 1 and Figure 2(e); note the difference in color scale between different velocity plots in Figure 2). To obtain the same jet velocity as in the nonmagnetic cases, much stronger forcing would be required. By implication, the amount of energy dissipation would be much higher. With a magnetic field amplitude B_o = 1 (case 3), which is equivalent to an Elsasser number Λ = 1 being reached for the conductivity at the bottom boundary, the jet velocity hardly changes with depth all the way down to the lower boundary. The Elsasser number is a measure for the ratio between Lorentz forces and Coriolis forces, and is defined as

$$\begin{eqnarray}&&{\rm{\Lambda }}=\displaystyle \frac{\sigma {B}^{2}}{\rho {\rm{\Omega }}}.\end{eqnarray} \tag{ 15 }$$

When the magnetic field strength is increased to B_o = 5 (case 4), with Λ = 25 being reached at the bottom, the surface wind speed inside the tangent cylinder drops by another factor of four compared to the case with B_o = 1. In contrast to the case with the weaker field, with B_o = 5 the wind velocity changes with depth in the lowermost part of the shell and drops by an order of magnitude toward the bottom (Figure 3). The same effect would be achieved for B_o = 1 by increasing the conductivity by a factor of 25. Note that values of Λ > 1 are not reached in the semiconducting regions in the gas planets, but require fully metallic conductivity. The zonal flow induces a magnetic field in the ϕ-direction (toroidal field, Figure 2(f)) in the lowermost part of the shell. However, the currents associated with this field, which run in the (s, z) plane, are not directly related to the decrease of the zonal velocity, which requires a j_ϕ current (Equation (14)). Instead, the Lorentz force associated with the toroidal field currents balances the Coriolis force of the meridional circulation in this depth range according to Equation (13). Figure 3 shows for the center of one of the jets inside the tangent cylinder, at s = 3, the variation of the zonal velocity and of the s component of the velocity as function of the z coordinate. Also shown is the driving force F_D and the ϕ component of the Lorentz force, taken from the simulation. Both forces are multiplied by 2/E_κ in the plot. They fall almost exactly on the curve for u_s, showing that the viscous term in Equation (13) is insignificant in this simulation. Compared to the previously discussed cases, where viscous forces balanced the meridional return flow, this role is now taken by magnetic forces. Finally, the interaction of the meridional flow with the magnetic field then gives rise to an electrical current j_ϕ, which relates to the drop of the zonal flow velocity with depth according to Equation (14). We also note that the Elsasser number, which varies strongly in the shell (mainly because of the variable conductivity), has a value around one at the point of the strongest decrease of u_ϕ(z).

Zoom In Zoom Out Reset image size

To achieve the same zonal flow velocity inside the tangent cylinder as in the nonconducting cases 1 and 2, the driving force—and hence, the power input and the dissipation—had to be increased by factors of approximately 100 for B_o = 1 and 500 for B_o = 5. At least in the quasi-linear regime of these simulations, where the imposed background field is not significantly altered, strong zonal flows do not abate drastically before they reach a depth of high conductivity. This seems incompatible with the secular variation and possibly the magnetic field geometry of Jupiter. An additional mechanism is needed to reduce the wind speed with depth.

While the simulations discussed so far confirm and illustrate some principles that had been conjectured before and studied in other models (e.g., Showman et al. 2006; Liu et al. 2008), we now enter new terrain by combining a stable layer with magnetic effects. Case 5, where B_o = 1, is shown in Figures 2(g)–(h). The zonal flow inside the tangent cylinder has a slightly larger surface velocity than in the nonmagnetic case without stable layer. In contrast to that case, the zonal flow velocity in case 5 starts to drop at mid-depth and decreases most sharply in the upper part of the stable layer (Figure 4(a)). Its value close to the inner boundary is approximately 1/1000th of the surface velocity. This drop-off is exclusively related to a thermal wind effect, as shown in Figure 4(b), which compares $\partial {u}_{\phi }/\partial z$ with the prediction from a thermal wind balance, i.e., with RaPrE_κ/(2r) ∂c/∂θ at s = 3. The two lines plot almost exactly on each other, showing that neither the Lorentz force term in Equation (14), which is essential in case 4 without stable layer, nor the viscous term play a role. However, the Lorentz force is important for the meridional circulation, which causes the codensity perturbation that is essential for the drop-off with depth of the jet velocity. In Figure 4(c), we compare the velocity in the s-direction on the cord with s = 3 with the various force terms in Equation (13). The Lorentz force is not the only contributor at the depths where the driving force F_D has dropped to zero, but is essential to balance F_D in the meridional flow equation. Without it, in case 2, where viscosity is the only antagonist to the driving force, the meridional circulation is unable to penetrate deeply into the stable layer and to create a codensity anomaly that kills the jets before they can reach the bottom boundary. However, in the magnetic case with a stable layer, the viscous force in Equation (13) is of the same magnitude as the Lorentz force. This contrasts with the strongly magnetic but neutrally stratified case 4, where the viscous contribution is insignificant. This difference is not surprising, as u_ϕ is almost three orders of magnitude smaller in the neutrally stratified case. That viscous effects are rather large is also shown by the fact that dissipation by ohmic losses is generally smaller than viscous dissipation, especially also in case 5 (D_ohm in Table 1). In a real planet, the opposite might be expected.

Zoom In Zoom Out Reset image size

All cases discussed so far have been calculated in an almost linear regime. The Rossby number $\mathrm{Ro}={U}_{\phi ,\max }/({\rm{\Omega }}D)$, which describes the ratio between the nonlinear inertial term and the Coriolis term, reaches only up to 2 × 10⁻⁴. In the magnetic cases, the amplitude of the induced field is much smaller than that of the imposed field, and the Lorentz force arising from the interaction of the induced currents with the induced field can be ignored. The Rossby number for Jupiter's zonal jets at higher latitudes is (1–2) × 10⁻². This Rossby number is reached in case 6, where the forcing amplitude A_f has been increased by a factor of 100 compared to that in case 5. Here, the amplitude of the induced toroidal field reaches 17% of that of the imposed dipole field. Nonetheless, all properties of the solution still scale almost linearly compared to case 5, i.e., they are larger by a factor of 100 (Table 1). The depth and the shape of the drop-off of the jet flow are essentially the same as in case 5. This suggests that linearized models, such as those studied in Section 4, can be meaningful for the gas planets.

To summarize the results of this section, only the combination of a stably stratified layer with electromagnetic effects seems to be able to simultaneously meet two requirements: (1) for a given moderate forcing, strong zonal jets develop at the surface; and (2) at the same time, their velocity drops off to very small values before the jet has reached a depth with high electrical conductivity. While the axisymmetric models are useful to illustrate the main physical effects, their excessive viscosity prevents them from being used directly for a quantitative comparison with the gas planets. This would require much lower values of the Ekman number that cannot be reached, for technical reasons, in such simulations.

The model setup is illustrated in Figure 5. A layer of unit depth (0 < z < 1) is assumed. Gravity and rotation are aligned and parallel to the z-axis. The imposed magnetic field ${{\boldsymbol{B}}}_{{\boldsymbol{o}}}$ is constant and points in the z-direction. This would approximately represent the situation close to the poles in a spherical shell. The results of the axisymmetric shell models suggest that the behavior is not fundamentally different at lower latitudes inside the tangent cylinder. Cartesian models with the magnetic field pointing in the x-direction (not discussed in detail here) show similar results. The system is periodic in the x-direction, which replaces latitude or the s-coordinate in a sphere. A flow in the y-direction (which replaces longitude) is driven by an imposed sinusoidal force F_D = A_f f(z) sin(kx). The system is two-dimensional, i.e., ∂/∂y = 0. Here, f(z) has the same form as in Equation (10), with z replacing r. Forcing is nonzero in the range 0.7 < z < 1.0.

Zoom In Zoom Out Reset image size

Equations (1)–(4) are written in terms of the following dependent variables: U = u_y for the jet flow, the stream function ψ for the meridional flow that is obtained as $\tilde{\rho }{{\boldsymbol{u}}}_{\mathrm{meri}}={\boldsymbol{\nabla }}\times (\tilde{\rho }\,\psi \,{{\boldsymbol{e}}}_{{\boldsymbol{y}}})$, c for codensity perturbation, b_z for the perturbation of the z-component of the poloidal field, and j_z for the z-component of the current density that is associated with the toroidal field. In contrast to the Boussinesq spherical shell models, we use here the anelastic approximation. For simplicity, we assume an exponential density variation with a constant scale height d_ρ.

The equations are linearized, assuming for the perturbation field $| {\boldsymbol{b}}| \ll {B}_{o}$ and for the meridional flow $| {u}_{x}| ,| {u}_{z}| \ll | U|$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial U}{\partial t}-\displaystyle \frac{2}{{E}_{\kappa }}\left(\displaystyle \frac{\partial \psi }{\partial z}-{d}_{\rho }^{-1}\psi \right)\\ \ \ =\Pr \left({{\boldsymbol{\nabla }}}^{2}U-{d}_{\rho }^{-1}\displaystyle \frac{\partial U}{\partial z}\right)+\displaystyle \frac{1}{{{qE}}_{\kappa }\tilde{\rho }}{B}_{o}\displaystyle \frac{\partial {b}_{y}}{\partial z}+{F}_{D},\end{array}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{\partial ({{\rm{\nabla }}}^{2}\psi -{d}_{\rho }^{-1}\partial \psi /\partial z)}{\partial t}-\displaystyle \frac{2}{{E}_{\kappa }}\displaystyle \frac{\partial U}{\partial z}\\ \ \ =\Pr \,{F}_{\nu }^{\psi }+\displaystyle \frac{1}{{{qE}}_{\kappa }\tilde{\rho }}{B}_{o}({{\rm{\nabla }}}^{2}{b}_{x}+{d}_{\rho }^{-1}{j}_{y})-\mathrm{RaPr}\displaystyle \frac{\partial c}{\partial x},\end{array}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial c}{\partial t}={{\rm{\nabla }}}^{2}\,c-{d}_{\rho }^{-1}\displaystyle \frac{\partial c}{\partial z}-\displaystyle \frac{\partial \psi }{\partial x}{C}_{o}^{{\prime} },\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {b}_{z}}{\partial t}-{B}_{o}\displaystyle \frac{{\partial }^{2}\psi }{\partial x\partial z}={q}^{-1}\eta {{\rm{\nabla }}}^{2}{b}_{z},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {j}_{z}}{\partial t}-{B}_{o}\displaystyle \frac{{\partial }^{2}U}{\partial x\partial z}={q}^{-1}\left(\eta {{\rm{\nabla }}}^{2}{j}_{z}+\displaystyle \frac{\partial \eta }{\partial z}\displaystyle \frac{\partial {j}_{z}}{\partial z}\right).\end{eqnarray} \tag{ 20 }$$

The viscous term in Equation (17) is

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\nu }^{\psi } & = & -{{\rm{\nabla }}}^{4}\psi +3{d}_{\rho }^{-1}{{\rm{\nabla }}}^{2}\displaystyle \frac{\partial \psi }{\partial z}-{d}_{\rho }^{-2}\left(3\displaystyle \frac{{\partial }^{2}\psi }{\partial {z}^{2}}+\displaystyle \frac{{\partial }^{2}\psi }{\partial {x}^{2}}\right)\\ & & +\,{d}_{\rho }^{-3}\displaystyle \frac{\partial \psi }{\partial z}.\end{array}\end{eqnarray} \tag{ 21 }$$

The components b_x, b_y that appear in these equations are obtained from the primary variables b_z, j_z by using ${\boldsymbol{\nabla }}\cdot {\boldsymbol{b}}=0$, ${\boldsymbol{\nabla }}\times {\boldsymbol{b}}={\boldsymbol{j}}$, and ${\boldsymbol{\nabla }}\cdot {\boldsymbol{j}}=0$.

The variables are expanded in sine or cosine functions in the x-direction, e.g., $U=\bar{U}(z)\sin ({kx})$. Because of the linearization, all wavenumbers k decouple. Equations (16)–(20) reduce to a set of five coupled ordinary differential equations in z. The same boundary conditions as in the spherical models are assumed, except for the magnetic field on the lower boundary, which is matched with a field that is obtained for an exponentially decreasing diffusivity η in the substratum z < 0. This condition can be implemented analytically for a given horizontal wavenumber k. The background codensity gradient changes around z = z_s from stable (${C}_{o}^{{\prime} }=1$) in the lower part to neutral above. The transition is smoothed again, but the transition interval is narrower than in the spherical shell simulations and restricted to z_s ± 0.01. The reference value for codensity diffusivity κ_o is now that in the lower layer. We assume that, in the upper region, the effective diffusivity by turbulent convection is so efficient that it practically wipes out any significant codensity differences. This means codensity is not solved for in that region, and the condition c = 0 is applied at the upper end of the stable region at z = z_s + 0.01.

We note that, in the case of stationary solutions (∂/∂t = 0), which are the only ones considered here, the Roberts number q can be eliminated from the equations by introducing a modified variable for the magnetic field perturbation and the current density, ${{\boldsymbol{b}}}^{* }={\boldsymbol{b}}/q$ and ${{\boldsymbol{j}}}^{* }={\boldsymbol{j}}/q$. For the simulations, where we set q = 1, this is not important. However, it becomes relevant when we scale the model results to real planets.

Equations (16)–(20) are solved for a specific wavenumber k, where derivatives ∂/∂x turn into factors ik, by a finite-difference scheme, using up to 4000 equidistant points in z. The equations are marched in time using an implicit scheme until a steady state has been reached. The solutions are linear in A_f, and therefore the forcing amplitude can be easily adjusted such that all reported solutions have a velocity of U = 1.0 at z = 1. We report results for E_κ in the range between 10⁻⁸ and 10⁻¹¹, and Rayleigh numbers Ra between 10¹⁶ and 3 × 10²⁴. In the cases with the highest values of the Rayleigh number, i.e., Ra ≥ 10²⁴ for z_s = 0.7 and Ra ≥ 10²¹ at z_s = 0.6, the numerical scheme becomes unstable for our standard setup. A stable solution can be obtained by applying the boundary conditions for U, ψ, and c not at z = 0 but at z = 0.35, and by setting these variables to zero below that depth. This is justified because at the highest values of Ra, all three variables assume almost zero values in the lower part of the domain. Comparison to a case in which a converged solution could also be obtained for the standard setup showed only very slight differences.

We strive at using parameters that are close to conditions in Jupiter as far as possible. We assume that the model domain covers the outer 10% of Jupiter, i.e., z = 0 corresponds to 0.9r_J. We use a superexponential variation of electrical conductivity. Figure 6 shows the conductivity values in the upper 12% of Jupiter's interior according to quantum mechanical ab initio simulations by French et al. (2012). They are decently fitted by

$$\begin{eqnarray}&&\sigma (r)=1.2\times {10}^{12}\,{\rm{S}}\,{{\rm{m}}}^{-1}\exp \left(\displaystyle \frac{-2.5}{1.05-r/{r}_{{\rm{J}}}}\right).\end{eqnarray} \tag{ 22 }$$

Identifying a radius of 0.9r_J with the lower boundary of the Cartesian model, Equation (22) translates into a nondimensional magnetic diffusivity of

$$\begin{eqnarray}&&\eta (z)=5.557\times {10}^{-8}\exp \left(\displaystyle \frac{25}{1.5-z}\right),\end{eqnarray} \tag{ 23 }$$

which is shown by the blue line in Figure 1(b) and is applied in the Cartesian models. To avoid numerical problems, the values are capped at a diffusivity ${\eta }_{\max }={10}^{10}$. Equation (22) gives a conductivity of 6.9 × 10⁴ S m⁻¹ at 0.9r_J. Together with a density at this depth of 577 kg m⁻³ (Nettelmann et al. 2012) and Ω = 1.76 × 10⁻⁴ s⁻¹, this results in a magnetic field scale of 1.2 mT. Jupiter's dipole field has, on average, a strength of 0.6 mT in the region (0.9–1.0)r_J. Hence, a nondimensional value of B_o = 0.5 is appropriate and is applied here.

Zoom In Zoom Out Reset image size

The average distance between a prograde jet and its retrograde neighbors at the surface of Jupiter is approximately 4° or 4800 km. With the depth of the simulation box corresponding to 0.1r_J or 7000 km, the nondimensional wavenumber is around $k=1.5\pi$, and we fix this value in our calculations.

We set the nondimensional density $\tilde{\rho }$ equal to one at z = 0. The scale height is quite variable in the molecular hydrogen envelope of the gas planets. For our calculations, where a constant scale height is assumed, we pick a nondimensional value of d_ρ = 0.43429, corresponding to a density ratio Δρ from top of bottom of ten. This roughly matches the scale height in Jupiter at about (0.95–0.96)r_J, or z = 0.5–0.6 in the simulations box, which is not far from where the drop-off of the jet flow occurs in our models.

We compare results for two locations of the upper boundary of the stable layer, z_s = 0.6 and z_s = 0.7, corresponding to 0.96r_J and 0.97r_J, respectively. As we will see, a stable region that extends rather high up turns out to be essential for meeting constraints on the power that is driving zonal winds in Jupiter. We calculated a series of models for each value of z_s by varying the Rayleigh number over a wide range, and choosing the Ekman number E_κ such that the fraction of viscous dissipation of the total (viscous plus ohmic) dissipation is always below 5%. At higher Rayleigh number, this requires lower values of the Ekman number.

Figure 7 shows depth profiles for a case with z_s = 0.7 at a high value of the Rayleigh number and low Ekman number. The jet velocity (full line in Figure 7(a)) stays constant above the stable region and drops sharply at the very top of the stable layer to very small values. The meridional circulation (u_x, u_z) is approximately seven orders of magnitude weaker than the jet flow. It vanishes at depths not far below the top of the stable layer. As was already seen in the spherical shell simulations, the decrease of U = u_y is associated with the thermal wind term F_B (the last term in Equation (17)), whereas the magnetic term F_L does not contribute (Figure 7(b)). The codensity perturbation causing the thermal wind effect is created by the meridional flow, which results from the Lorentz force associated with the induced electrical current component j_x as the only significant antagonist to the driving force F_D (Figure 7(c)). This contrasts with the spherical shell model, where viscous friction still played a significant role in the meridional force balance (compare Figure 4(c)).

Zoom In Zoom Out Reset image size

The electrical current (j_x, j_z) associated with the toroidal magnetic field is induced by the interaction of the zonal flow U with the background field in the region where U has not yet vanished, and where at the same time the electrical conductivity is already sufficient to allow for a current to develop (Figure 7(d)). This region is restricted to a rather thin sublayer near the top of the stable region. The electrical current closes via the j_z component through the highly conducting interior. Because of the high value of the Rayleigh number, a very weak codensity perturbation is sufficient to balance the drop of the zonal flow velocity via the thermal wind term. To set up the codensity perturbation, a weak meridional flow is sufficient. This, in turn, requires only a weak Lorentz force that is achieved by the small current j_x, which is possible in a region of poor conductivity. The Elsasser number Λ, which for a given B_o may be considered as a measure for the conductivity, is only 7 × 10⁻⁷ at the depth z_s = 0.69 where the sharpest drop in velocity occurs. This is in drastic contrast to the value Λ ≈ 1 that must be reached to let the zonal velocity drop in the absence of stable stratification (see Section 3.3). The ohmic dissipation associated with the current system is also concentrated around the transition between the stable and the unstable layer (red line in Figure 7(d)). Because of the weakness of the current, the dissipation remains comparatively small despite the high resistivity. Correspondingly, only a small driving power is required.

While the current components j_x, j_z associated with the induced toroidal magnetic field are already small, the j_y component that relates to the perturbation of the poloidal magnetic field is even smaller by several orders of magnitude (broken line in Figure 7(d)). Consequently, the b_x and b_z components are very tiny and the Lorentz force term associated with b_x in Equation (17) can be neglected.

Next, we want to quantify the dependence of several properties of interest on the control parameters. One such property is the power exerted by the driving forces, or equivalently, the energy dissipation. The nondimensional power per unit area of the surface, scaled by ${\rho }_{o}{\kappa }_{o}^{3}/{D}^{3}$, is obtained as

$$\begin{eqnarray}&&P=\displaystyle \frac{1}{2}{\int }_{0}^{1}\tilde{\rho }{F}_{D}(z)U(z){dz}.\end{eqnarray} \tag{ 24 }$$

The prefactor comes from averaging the sinusoidal variation in x. Given that U ≈ 1 in the depth range where F_D does not vanish, for the chosen forcing function and density profile, the power is in close approximation P = 0.00935A_f. Another property of interest is the characteristic length scale on which the jet velocity decreases. It is defined here as

$$\begin{eqnarray}&&L=\displaystyle \frac{5}{4}({z}_{90}-{z}_{10}),\end{eqnarray} \tag{ 25 }$$

with z₉₀ and z₁₀ being the points where U has dropped to 90% and 10%, respectively, of its surface value. A different definition of the decay length based on the value of dU/dz at the depth where U has dropped to 1/2 gives results that are not much different, but systematically smaller by a factor of about 0.85. The final property whose systematic dependence on control parameters we consider is the maximum strength of the toroidal magnetic field ${b}_{y,\max }^{* }$.

To identify the key parameter that controls these properties, we rewrite Equations (16)–(20) for a particular wavenumber k, drop the time derivative and the viscous terms, and eliminate the Roberts number q by using the modified variables for the perturbed field and for the current density. Furthermore, we drop the Lorentz force term in Equation (17), which was found to be insignificant. Because of the latter, we can ignore Equation (19), which is then only "diagnostic" and does not influence the other variables:

$$\begin{eqnarray}&&-2\left(\displaystyle \frac{d\psi }{{dz}},-,{d}_{\rho }^{-1},\psi \right)=-i\displaystyle \frac{{B}_{o}}{k}\displaystyle \frac{{{dj}}_{z}^{* }}{{dz}}+{\bar{A}}_{f}f(z),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&2\displaystyle \frac{{dU}}{{dz}}={ik}\,{\mathrm{Ra}}_{{\rm{\Omega }}}\,c,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\left(\displaystyle \frac{{d}^{2}}{{{dz}}^{2}}-{d}_{\rho }^{-1}\displaystyle \frac{d}{{dz}}-{k}^{2}\right)c={ik}\psi {C}_{o}^{{\prime} },\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&-{{ikB}}_{o}\displaystyle \frac{{dU}}{{dz}}=\eta \left(\displaystyle \frac{{d}^{2}}{{{dz}}^{2}}-{k}^{2}\right){j}_{z}^{* }+\displaystyle \frac{d\eta }{{dz}}\displaystyle \frac{{{dj}}_{z}^{* }}{{dz}}.\end{eqnarray} \tag{ 29 }$$

Here, we set ${E}_{\kappa }{A}_{f}={\bar{A}}_{f}$, and we define a "rotational Rayleigh number"

$$\begin{eqnarray}&&{\mathrm{Ra}}_{{\rm{\Omega }}}=\mathrm{Ra}{E}_{\kappa }\Pr ={\tilde{N}}^{2}/{E}_{\kappa }=\displaystyle \frac{{{gC}}_{o}^{{\prime} }{D}^{2}}{{\kappa }_{o}{\rm{\Omega }}}.\end{eqnarray} \tag{ 30 }$$

For fixed values of B_o and k, and fixed η(z), ${C}_{o}^{{\prime} }(z)$, $\tilde{\rho }(z)$ and forcing function f(z), Ra_Ω is the only remaining free parameter. Because we require U to become one at the surface, the forcing amplitude A_f or ${\bar{A}}_{f}$ is not a free control parameter but a diagnostic parameter, i.e., it sets the force needed to obtain a jet with unit velocity. Because the power P is proportional to A_f, it follows that the product E_κP is a function of Ra_Ω. Furthermore, the other diagnostic parameters, L and ${b}_{y,\max }^{* }$ (=${j}_{z,\max }^{* }/k$), are only functions of Ra_Ω. For a few cases, this has been verified by changing E_κ while keeping the same value of Ra_Ω. Aside from tiny differences due to residual viscosity effects, the results are identical.

With increasing values of Ra_Ω, the power required for driving a zonal flow with unit velocity at the surface decreases significantly (Figure 8(a)). The drop-off of the jet velocity near the top of the stable layer becomes sharper (Figure 8(b)), and the induced toroidal magnetic field becomes weaker (Figure 8(c)). The variation is not linear in the log–log plot, i.e., no simple power-law dependence applies. The curves flatten out toward high values of Ra_Ω for power and toroidal field amplitude. This is probably a consequence of the fact that, with increasing values of Ra_Ω, the action is squeezed more and more into a thin layer near the top of the stable region (Figure 7). For L, the curve steepens toward the high end of parameter values.

Zoom In Zoom Out Reset image size

While the drop-off length scale L does not differ much between the two model series, the driving power and the toroidal field strength depend rather strongly on the location of the upper boundary of the stable region, especially near the high end of values for Ra_Ω. For z_s = 0.7, the toroidal field is weaker and the power requirement is lower by more than an order of magnitude than it is for z_s = 0.6.

In order to apply the results to the gas planets, we need to estimate the appropriate values of the control parameters E_κ, q, and ${\mathrm{Ra}}_{{\rm{\Omega }}}={\tilde{N}}^{2}/{E}_{\kappa }$. We do this for Jupiter. The relevant value of the diffusivity κ is an effective one in the stable region. The molecular values of thermal and compositional diffusivity are not very different from each other in the range (0.25–1) × 10⁻⁶ m² s⁻¹ (French et al. 2012). For the nominal depth D = 7000 km of the Cartesian models, this would lead to values of E_κ on the order of 10⁻¹⁶. However, a stable region in a gas planet is unlikely to be stagnant. Transport processes occur by double-diffusive or "semi"-convection (Leconte & Chabrier 2013; Debras & Chabrier 2019), e.g., in the form of a stack of internally convecting layers. If Jupiter's internal heat flow of 5.4 W m⁻² (Guillot et al. 2004) had to be transported through the stable layer by conduction with the molecular conductivity λ ≈ 2 W m⁻¹ K⁻¹, it would require an unreasonable thermal gradient of 2700 K km⁻¹, very much larger than the adiabatic temperature gradient of 0.7 K km⁻¹. If we assume that the mean temperature gradient across the stack of double-diffusive layers is 1 K km⁻¹, i.e., moderately superadiabatic, the "effective" diffusivity due to the small-scale convection in the layers would be enhanced above the molecular value by a factor of ~2500. Applying this factor, we assume a value of 1.5 × 10⁻³ m² s⁻¹ for the effective diffusivity κ_eff. The effective Ekman number E_κ is then 1.7 × 10⁻¹³. Together with a plausible range of $0.8\lt \tilde{N}\lt 3$ for a stable layer in Jupiter suggested by the models of Debras & Chabrier (2019), this leads to values around 10¹³ for the parameter Ra_Ω. Relevant dimensional and resulting nondimensional parameters are summarized in Table 2.

Table 2.  Jupiter Parameters

D	${\rho }_{o}$	ηo	κeff	Ω	N/Ω	Eκ	q	Pr	RaΩ
7 × 106	577	12	1.5 × 10−3	1.76 × 10−4	0.8–3	1.7 × 10−13	1.2 × 10−4	≈1	4 × 1012–5 × 1013

Note. SI units where applicable.

Download table as:  ASCIITypeset image

The most extreme cases in the Cartesian model series reach the plausible range of Ra_Ω in Jupiter, shown as the shaded region in Figure 8. If we assume that the stable layer extends out to 0.97r_J (z_s = 0.7), our results predict for Jupiter a nondimensional power E_κP around 2 × 10⁻⁹, a drop-off length scale L = 0.025–0.04, and a toroidal field ${b}_{y}^{* }={b}_{y}/q$ of the order (0.7–1) × 10⁻⁸. These values refer to a nondimensional jet velocity of one. For a typical jet velocity inside Jupiter's tangent cylinder of 20 m s⁻¹, the nondimensional value is 10¹¹, using the effective value of κ. The result for the toroidal field strength is to be multiplied with this factor and that for power with the factor squared. We then end up with the following estimates for the dimensional values: power per unit area is around 0.7 W m⁻², the toroidal field strength is 0.11–0.15 mT, and the decay length of the jet flow is 200–300 km.

It is reassuring that the estimated power that is pumped into the zonal jets is less than the internal heat flow by almost an order of magnitude. To convert such fraction of the heat flow into kinetic energy poses no problem from a thermodynamic efficiency point of view. The analysis of the flux of momentum from eddies into zonal flow, based on cloud-tracking data at Jupiter, suggests that the associated energy conversion rate could be 0.7–1.2 W m⁻² (Salyk et al. 2006), in full agreement with our result. The predicted toroidal field strength is one-fifth of that of the dynamo-generated dipole field. Therefore, the linearity assumption in the calculations is not violated—albeit, it is just marginally satisfied.

These estimates hold for a location of the upper boundary of the stable layer at 0.97r_J. Figure 8 shows that lowering the boundary to 0.96r_J would boost the power requirement by a factor of 30, which may be incompatible with the probable upper limit. This power estimate must be taken with caution, because also the maximum toroidal field gets larger by an order of magnitude in the simulations. When scaled to Jupiter parameters, it would exceed the imposed background field. This agrees with the finding of Wicht et al. (2019b), that the wind-induced zonal toroidal field reaches the same order of magnitude as the poloidal field for a penetration depth of the winds between 0.97r_J and 0.96r_J. Hence, the linearity assumption of our simulations would be violated. Nonetheless, it may be safe to take a value around 0.965r_J as the plausible lower limit for the top of the stable layer.