THE SUN'S MERIDIONAL CIRCULATION AND INTERIOR MAGNETIC FIELD

We measure distances in units of the solar radius, R_☉ ≃ 7 × 10¹⁰ cm, and velocities in units of R_☉Ω_☉, where Ω_☉ ≃ 3 × 10⁻⁶ s⁻¹ is the mean solar rotation rate, which we use as the inverse timescale. As illustrated in Figure 1, the Cartesian coordinate system (x, y, z) is chosen such that x is the azimuthal coordinate, y ∈ [0, π] is the latitudinal coordinate, and z ∈ [0, 1] is the vertical coordinate. Although there is no precise equivalence between our Cartesian geometry and the spherical geometry of the solar interior, for the sake of interpretation we will take y = 0 and y = π to be the "poles" and y = π/2 to be the "equator." The radiative–convective interface is at z = h ≃ 0.7, so that h < z < 1 represents the convection zone, and 0 < z < h represents the radiative zone.

Zoom In Zoom Out Reset image size

The steady background state considered is depicted in Figure 1. We model the transition between the convective and radiative zones by imposing a vertical profile for the dimensionless buoyancy frequency n(z) with

$$\begin{equation} n^2(z)=\frac{n^2_{\rm rz}}{2}\left\lbrace 1+\tanh \left(\frac{h-z}{\Delta }\right)\right\rbrace, \end{equation} \tag{ 4 }$$

so that n = 0 in the convection zone and n = n_rz > 0 in the radiative zone. The lengthscale Δ is introduced in order to ensure the smoothness of the background state, which is necessary for numerical reasons. We may regard Δ as the thickness of the convective overshoot region at the base of the convection zone.

We impose a horizontal background magnetic field B₀ of the form

$$\begin{equation} {\bf B}_0 = B_0(z)\hat{{\bf e}}_y, \end{equation} \tag{ 5 }$$

where $\hat{{\bf e}}_y$ is the unit vector in the y direction. In order to confine this field against ohmic diffusion, in the steady state, we must invoke a background meridional flow. The simplest case is that of uniform downwelling of dimensionless magnitude U₀, and

$$\begin{equation} B_0(z) \;=\; {\cal B}_0 {e}^{-z/\delta }, \end{equation} \tag{ 6 }$$

where ${\cal B}_0$ is a constant and δ is the dimensionless lengthscale

$$\begin{equation} \delta \;=\; \frac{\eta /\big(R_\odot ^2\Omega _\odot\big)}{U_0} = \frac{{E}_\eta }{U_0}. \end{equation} \tag{ 7 }$$

Here, η and E_η represent the magnetic diffusivity in dimensional and non-dimensional units, respectively; we call E_η the "magnetic Ekman number." Although this background state cannot hold over all latitudes and depths within the solar interior (for reasons of mass conservation and solenoidality of B₀), we hope that it offers a reasonable qualitative approximation to conditions within the tachocline at middle and high latitudes.

We note that B₀ is not a force-free field since it has a non-constant magnetic pressure B²₀/(8π). However, the magnetic pressure can be balanced by a perturbation to the background gas pressure p₀. Estimates of the strength of the Sun's interior magnetic field typically lie within the range 10⁻³ to 10³ G (e.g., Mestel & Weiss 1987; Gough & McIntyre 1998; see also Section 4.3). For field strengths in this range, the gas pressure far exceeds the magnetic pressure, and so the required perturbation is very small.

In what follows we consider linear, steady-state perturbations to this background state. The principle shortcoming of our model is, arguably, the artificial manner in which the magnetic field is confined to the interior via the imposed downwelling U₀. Magnetic confinement is thereby built into the model and controlled by the parameters ${\cal B}_0$ and δ. This idealization is necessary to allow a linear study. A complete model of the tachocline would need to describe self-consistently the processes that act to confine the magnetic field, rather than treating it as part of the background. Such a model would necessarily be nonlinear and awaits future numerical work. We note also that our Cartesian model cannot adequately describe the polar regions, where the effects of geometrical curvature become significant. The polar magnetic confinement problem has been studied in detail elsewhere (Wood & McIntyre 2011), and self-consistent, fully nonlinear solutions have been obtained.

As in the studies of GAA09 and GB10 we do not model the turbulence in the convection zone in detail. Instead, we introduce a simple parametric model to describe its role in driving a large-scale differential rotation profile and in enhancing the transport of heat. In the region z ∈ [h, 1], we model the turbulent heat transport as a diffusive process and replace the divergence of the Reynolds stress with a forcing term that represents linear relaxation toward a prescribed rotation profile, which we choose to mimic the observed differential rotation of the solar convection zone.

For simplicity, we also make the Boussinesq approximation (e.g., Spiegel & Veronis 1960), in which we neglect variations in the background pressure p₀, density ρ₀, and temperature T₀, and we neglect pressure perturbations in the equation of state. Although the Boussinesq approximation cannot be rigorously justified over the entirety of our domain, we do not expect this approximation to have a qualitative effect on our results. Indeed, we are interested primarily in solutions for which the meridional flows and differential rotation are confined to a thin, tachocline-like layer below the convection zone, within which the Boussinesq approximation should be quite accurate.

The non-dimensional governing equations, linearized about the background state described in the previous section, are

$$\begin{eqnarray} 2\hat{{\bf e}}_z\times {\bf u}&=& {-}\nabla p + T \hat{{\bf e}}_z + \Lambda {e}^{-z/\delta }\left[ \frac{1}{\delta }\hat{{\bf e}}_x \times {\bf B}+ (\nabla \times {\bf B}) \times \hat{{\bf e}}_y \right]\nonumber\\ && + {E}_\nu \nabla ^2{\bf u}- \lambda (z)({\bf u}-{\bf u}_{\rm cz}) \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} n^2(z)\, w &= {E}_\kappa \nabla \cdot (f(z)\,\nabla T) \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} 0 &= \nabla \times \left({\bf u}\times {e}^{-z/\delta }\hat{{\bf e}}_y - \frac{1}{\delta }\hat{{\bf e}}_z \times {\bf B}- \nabla \times {\bf B}\right) \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \nabla \cdot {\bf u}&= 0 \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \nabla \cdot {\bf B}&= 0\,, \end{eqnarray} \tag{ 12 }$$

where u = (u_x, u_y, u_z) is the velocity perturbation, B = (b_x, b_y, b_z) is the perturbation to the magnetic field, p is the pressure perturbation, and T is the temperature perturbation. These equations have been scaled using ρ₀R²_☉Ω²_☉ as the unit of pressure, T₀R_☉Ω²_☉/g as the unit of temperature, and ${\cal B}_0/{E}_\eta$ as the unit of magnetic field strength. The following dimensionless parameters have been introduced:

$$\begin{eqnarray} \Lambda =& \frac{{\cal B}_0^2 }{4\pi \eta \rho _0\Omega _\odot }\,,& \;\; \mbox{the Elsasser number at $z=0$;} \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} {E}_\nu =& \frac{\nu }{R_\odot ^2\Omega _\odot }\,,& \;\; \mbox{the Ekman number;}\qquad\qquad\quad \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} {E}_\kappa =& \frac{\kappa }{R_\odot ^2\Omega _\odot }\,,& \;\; \mbox{the thermal Ekman number in the}\nonumber\\ &&{\rm radiative\; interior.} \end{eqnarray} \tag{ 15 }$$

Note that the non-dimensionalized perturbation equations have no explicit dependence on the magnetic Ekman number E_η. However, because E_η has been used to scale the magnetic field, its value governs the magnitude of perturbations relative to the background field B₀.

The Elsasser number Λ is a measure of the relative magnitude of Lorentz forces compared to Coriolis forces. It is convenient to define also a local Elsasser number,

$$\begin{eqnarray} \Lambda _{\rm loc}(z) &=& \frac{B_0^2(z)}{4\pi \eta \rho _0\Omega _\odot } \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&\quad\quad = \Lambda {e}^{-2z/\delta }\,, \end{eqnarray} \tag{ 17 }$$

in order to measure the relative magnitudes of these forces at each altitude z.

In both the momentum equation (8) and the thermal energy equation (9), we have neglected terms describing advection by the background downwelling U₀. It can be shown that, in the results presented here, including these terms would produce only a small correction to the solutions. Neglecting these terms simplifies the analysis and reduces the dependence of the solutions on the artificially imposed flow U₀.

The final term in the linearized momentum equation (8) describes a linear relaxation toward a prescribed rotation profile u_cz. The dimensionless parameter λ, which determines the rate of relaxation, can be regarded as the inverse of the convective turnover timescale. This parameterization for the angular momentum transport by Reynolds stresses is similar to that used by Bretherton & Spiegel (1968), GAA09, and GB10. For simplicity, we assume that λ takes the constant value 1/τ_c in the convection zone and that the prescribed differential rotation u_cz is sinusoidal in latitude. Following GAA09, we adopt the profiles

$$\begin{equation} {\bf u}_{\rm cz} = u_{\rm cz}(z)\, {e}^{{i}ky}\hat{{\bf e}}_x, \end{equation} \tag{ 18 }$$

$$\begin{equation} \mbox{and}\quad \lambda (z) = \frac{1}{2\tau _c}\left\lbrace 1+\tanh \left(\frac{z-h}{\Delta }\right)\right\rbrace, \end{equation} \tag{ 19 }$$

where Δ is the same "overshoot" length as in Equation (4). In the numerical results presented later we set τ_c = 0.1 and k = 2, representing equatorial symmetry at the solar surface, and u_cz(z) is taken to be linear in z. We show in Appendix B that the flow below the convection zone is not very sensitive to the value of τ_c, or to the form of u_cz(z), and that the results are easily extended to any particular u_cz(z). This insensitivity to the choice of parameterization within the convection zone reflects the robustness of both the gyroscopic pumping mechanism and the burrowing tendency of meridional flows.

Finally, we model turbulent heat transport within the convection zone as an increase in the thermal diffusivity by a factor f relative to its laminar, microscopic value. We take

$$\begin{equation} f(z) = 1 + \frac{f_0-1}{2}\left\lbrace 1+\tanh \left(\frac{z-h}{\Delta }\right)\right\rbrace \,, \end{equation} \tag{ 20 }$$

so that f  =  1 within the radiative zone and f  =  f₀  ≫  1 within the convection zone. Thermal perturbations at the radiative–convective interface are therefore communicated more efficiently into the convection zone than into the radiative zone. In the limit f₀ → ∞, we expect the convection zone to become isothermal, which is equivalent to isentropic in our Boussinesq model. In the numerical results presented in Section 4, we take f₀ = 10⁶. The physical significance of f₀, and its effect on the solutions, is described in Section 4.3, below Equation (72).

We assume that the viscosity ν and the magnetic diffusivity η retain their laminar, microscopic values throughout the solar interior. Introducing larger "turbulent" values for ν and η within the convection zone would not significantly affect the results we present here.

The perturbation fields must have the same sinusoidal dependence on y as the forcing term λ(z) u_cz in Equation (8), and so we seek solutions of the form

$$\begin{equation} q(y,z)=\hat{q}(z)\,{e}^{iky}\,. \end{equation} \tag{ 21 }$$

The perturbation equations (8)–(12) then become a system of ordinary differential equations for the perturbation amplitudes $\hat{q}$:

$$\begin{eqnarray} {-}2\hat{u}_y&=& {i}k \Lambda \hat{b}_x{e}^{-z/\delta }+{E}_\nu \left(\frac{{d}^2\hat{u}_x}{{d}z^2}-k^2\hat{u}_x\right) - \lambda (z)(\hat{u}_x-u_{\rm cz})\nonumber\\ \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} 2\hat{u}_x&=& {-}{i}k\hat{p} - \frac{\Lambda }{\delta }\hat{b}_z{e}^{-z/\delta } + {E}_\nu \left(\frac{{d}^2\hat{u}_y}{{d}z^2}-k^2\hat{u}_y\right) - \lambda (z)\hat{u}_y\nonumber\\ \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} \frac{{d}\hat{p}}{{d}z} - \hat{T} &=& \frac{\Lambda }{\delta }\hat{b}_y{e}^{-z/\delta } + \Lambda \left({i}k\hat{b}_z-\frac{{d}\hat{b}_y}{{d}z}\right){e}^{-z/\delta }\nonumber\\ &&+{E}_\nu \left(\frac{{d}^2\hat{u}_z}{{d}z^2}-k^2\hat{u}_z\right) - \lambda (z)\hat{u}_z \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} 0 &=& {i}k\hat{u}_x{e}^{-z/\delta } + \frac{1}{\delta } \frac{{d}\hat{b}_x}{{d}z} + \frac{{d}^2\hat{b}_x}{{d}z^2}-k^2\hat{b}_x \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} 0 &=& {i}k\hat{u}_z{e}^{-z/\delta } + \frac{1}{\delta } \frac{{d}\hat{b}_z}{{d}z} + \frac{{d}^2\hat{b}_z}{{d}z^2}-k^2\hat{b}_z \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} \frac{n^2(z)}{{E}_\kappa }\hat{u}_z&=& \frac{{d}}{{d}z}\left(f(z)\frac{{d}\hat{T}}{{d}z}\right) - f(z) k^2\hat{T} \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} {i}k\hat{u}_y+ \frac{{d}\hat{u}_z}{{d}z} &=& 0 \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} {i}k\hat{b}_y+ \frac{{d}\hat{b}_z}{{d}z} &=& 0 \,. \end{eqnarray} \tag{ 29 }$$

The system of Equations (22)–(29) is solved numerically using a two-point boundary-value solver, given the boundary conditions listed below.

The upper and lower boundaries are modeled as impenetrable, (viscous) stress-free, isothermal, and electrically insulating. That is, we impose

$$\begin{eqnarray} \left. \begin{array}{@{}rcl@{}}\hat{u}_z&=& 0 \\ \ms \dsty\frac{{d}\hat{u}_x}{{d}z} &=& 0 \\ \ms \dsty\frac{{d}\hat{u}_y}{{d}z} &=& 0 \\ \ms \hat{T} &=& 0 \\ \ms \hat{b}_x&=& 0 \end{array} \right\rbrace \;\;\; &\mbox{at} \;\;\; z = 0 \;\; \mbox{and} \;\; z = 1,& \end{eqnarray} \tag{ 30 }$$

$$\begin{equation} \frac{{d}\hat{b}_z}{{d}z} = k\hat{b}_z\;\;\; \mbox{at} \;\;\; z = 0, \end{equation} \tag{ 31 }$$

$$\begin{equation} \mbox{and}\quad \frac{{d}\hat{b}_z}{{d}z} = -k\hat{b}_z\;\;\; \mbox{at} \;\;\; z = 1. \end{equation} \tag{ 32 }$$

Most of these boundary conditions have been selected for simplicity and convenience. In fact, the solutions turn out to be rather insensitive to the choice of boundary conditions in the parameter regime of interest, as we have verified by comparing solutions obtained with different sets of boundary conditions.

In order to delineate the various physical effects present in the solutions, we begin by considering the unstratified case, n_rz = 0. Effects associated with stratification will be described in detail in Section 4.

A series of results is shown in Figure 2. In the absence of any magnetic field (i.e., setting Λ = 0), we recover the non-magnetic, unstratified solution of GAA09 (see their Figure 2). Since the lower boundary is stress-free, the differential rotation imposed within the convection zone is able to spread throughout the radiative zone and establishes a state close to Taylor–Proudman balance. In this steady state, the only meridional flows are those maintained by viscosity, and their magnitude scales with the Ekman number, E_ν. The meridional flows therefore become vanishingly weak in the limit E_ν → 0. When Λ > 0, on the other hand, the Lorentz force produces departures from Taylor–Proudman balance and maintains a meridional flow even in the limit E_ν → 0 (Figure 2(b)).

Zoom In Zoom Out Reset image size

For Λ of order unity, the Lorentz forces are most significant close to the lower boundary z = 0, within a boundary-layer analogous to an Ekman–Hartmann layer. The structure of the boundary layer depends on the particular choice of boundary conditions at z = 0, including the magnetic boundary conditions. However, we have verified that the meridional flow outside the boundary layer is hardly affected by the choice of magnetic boundary conditions.

For Λ ≫ 1 the radiative zone divides into two dynamically distinct regions, separated by a thin internal boundary layer at z = z₀, say (see Figure 2(c)). The region above the boundary layer is close to Taylor–Proudman balance, with all components of u independent of z, indicating that Lorentz forces are locally negligible (see Figure 8(a) in Section 4.1 for plots of $\hat{u}_x$). In the region below the boundary layer, on the other hand, the Lorentz force is dominant, and all components of u decay exponentially with depth. From here on we refer to these two regions as being "weakly magnetic" and "magnetically dominated," respectively. The boundary layer separating these two regions we call the "magnetic transition layer." Since the flow within the magnetically dominated region decays exponentially with depth, the solution is insensitive to the choice of boundary conditions at z = 0.

The structure of the weakly magnetic and magnetically dominated regions is illustrated further in Figure 3, which shows the meridional streamlines and magnetic field lines, as well as contours of the azimuthal components of the velocity and magnetic fields, in a vertical cross-section through a typical solution. Within the convection zone, the forcing that drives differential rotation also gyroscopically pumps a meridional flow. Below the base of the convection zone the forcing is switched off, and Taylor–Proudman balance is achieved. As a result, the meridional streamlines and angular velocity contours become aligned with the rotation axis, and a portion of the convection zone's differential rotation extends into the radiative zone. This differential rotation winds up the magnetic field lines, producing a Lorentz torque whose amplitude increases exponentially with depth. Within the magnetic transition layer (here around z₀ = 0.2), the Lorentz torque is strong enough to overcome Taylor–Proudman balance and gyroscopically pumps a meridional flow. The magnetic field lines are dragged upward where the flow is upwelling, 0  <  y  <  π/2, and pushed downward where the flow is downwelling, π/2  <  y  <  π (not shown in the figure). Below the transition layer, meridional flows and differential rotation are both strongly suppressed, and consequently perturbations to the magnetic field lines are smaller.

Zoom In Zoom Out Reset image size

Figure 4 shows the result of increasing Λ still further. Over many orders of magnitude the only effect of increasing Λ is to raise the position of the magnetic transition layer, shrinking the vertical extent of the weakly magnetic region. The dynamics within the convection zone remain unaffected until the point at which the transition layer meets the base of the convection zone.

Zoom In Zoom Out Reset image size

It is obviously of interest to identify the factors determining the location of the magnetic transition layer. We anticipate that this transition occurs where the local Elsasser number Λ_loc, defined in Equation (17), takes a critical value. We therefore expect the vertical position of the transition layer, z₀, to follow a law of the form

$$\begin{equation} z_0 = \mbox{const.} + \frac{1}{2}\delta \ln {\Lambda } \end{equation} \tag{ 33 }$$

as Λ is increased, where the constant depends on the aforementioned critical value of Λ_loc. It is readily verified from Figure 4 that z₀ does indeed follow such a law. In the following sections, we analyze the structure of the solutions in detail in order to determine both the critical value of Λ_loc and the amplitude of the gyroscopically pumped meridional flow below the convection zone.

In order to understand the results presented in Section 3.1 more quantitatively, we now derive approximate analytical solutions of the governing equations. In Section 3.2.1, we derive a boundary-layer solution for the magnetic transition layer. In Section 3.2.2, we use that boundary-layer solution to construct a global analytical solution.

3.2.1. Transition Layer Solution

We seek to describe the transition layer between the weakly magnetic and magnetically dominated regions of the flow. Since this transition occurs below the base of the convection zone we may neglect the terms in the momentum equation involving λ. Based on the numerical solutions shown above, we anticipate that the thickness of the transition layer is of the same order as δ, the scale height of the background magnetic field. Since δ ≪ 1/k, we make the boundary-layer approximation

$$\begin{equation} \left(\frac{{d}^2}{{d}z^2} - k^2\right) \simeq \frac{{d}^2}{{d}z^2}\,. \end{equation} \tag{ 34 }$$

We also neglect viscosity, i.e., we set E_ν = 0. It can be verified a posteriori that viscosity plays no significant role in the dynamics of the transition layer provided that δ is much larger than the Ekman length, that is, provided that

$$\begin{equation} \delta \gg {E}_\nu ^{1/2}. \end{equation} \tag{ 35 }$$

This condition is easily satisfied in the solar interior, where E_ν ≃ 10⁻¹⁵.

With these approximations, the governing equations (22)–(29) reduce to the following set:

$$\begin{eqnarray} -2\hat{u}_y&=& {i}k \Lambda \hat{b}_x{e}^{-z/\delta } \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} 2\hat{u}_x&=& -{i}k\hat{p}-\frac{\Lambda }{\delta } \hat{b}_z{e}^{-z/\delta } \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} \frac{{d}\hat{p}}{{d}z} &=& \frac{\Lambda }{\delta }\hat{b}_y{e}^{-z/\delta } - \Lambda \frac{{d}\hat{b}_y}{{d}z}{e}^{-z/\delta } \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} 0 &=& {i}k\hat{u}_x{e}^{-z/\delta } + \frac{1}{\delta }\frac{{d}\hat{b}_x}{{d}z} + \frac{{d}^2\hat{b}_x}{{d}z^2} \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} 0 &=& {i}k\hat{u}_z{e}^{-z/\delta } + \frac{1}{\delta }\frac{{d}\hat{b}_z}{{d}z} + \frac{{d}^2\hat{b}_z}{{d}z^2} \end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray} {i}k\hat{u}_y+\frac{{d}\hat{u}_z}{{d}z} &=& 0 \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} {i}k\hat{b}_y+\frac{{d}\hat{b}_z}{{d}z} &=& 0 \,. \end{eqnarray} \tag{ 42 }$$

We show in Appendix A that these can be combined into a single equation for $\hat{u}_z$:

$$\begin{equation} \left(\frac{{d}}{{d}z}-\frac{2}{\delta }\right) \left[ \left({e}^{z/\delta }\frac{{d}}{{d}z}\right)^4 + \frac{1}{4}k^4\Lambda ^2 \right]\hat{u}_z= 0\,. \end{equation} \tag{ 43 }$$

From inspection of Equation (43), we deduce that the vertical position of the transition layer, z₀, is given by

$$\begin{eqnarray} ({e}^{z_0/\delta }/\delta)^4 &= \frac{1}{4}k^4\Lambda ^2 \nonumber \\ \Rightarrow z_0 &= \delta \ln (k\delta /\sqrt{2}) + \frac{1}{2}\delta \ln {\Lambda }, \end{eqnarray} \tag{ 44 }$$

which is consistent with Equation (33). Furthermore, this expression shows that the transition layer is located where the local Elsasser number Λ_loc, defined by Equation (17), takes the critical value 2/(kδ)². This corresponds to a critical magnetic field strength

$$\begin{equation} B_{\rm crit} = \frac{\sqrt{8\pi \eta \rho _0\Omega _\odot }}{k\delta } \simeq 40 \left(\frac{\eta }{\eta _\odot } \right)^{1/2} \left(\frac{0.001}{\delta } \right) {\rm G} \,, \end{equation} \tag{ 45 }$$

where ρ₀ ≃ 0.2 g cm⁻³ and η_☉ ≃ 400 cm² s⁻¹.

We seek the solution to Equation (43) that has the correct behavior in both the weakly magnetic region, for z − z₀ ≫ δ, and the magnetically dominated region, for z − z₀ ≪ −δ. In particular, the solution must match onto the differential rotation, $\hat{u}_x= u_{\rm t}$ say, in the region immediately above the transition layer. Here, and subsequently, we use a subscript "t" to refer to the region immediately above the transition layer, which is also the bottom of the weakly magnetic region. We show in Appendix A that the unique solution with the required properties is

$$\begin{equation} \hat{u}_z= {i}k\delta u_{\rm t}\left[\raisebox {0.5cm}{} I_1(\zeta) - \frac{\pi }{2}\,{\rm Re}\left\lbrace \exp \left(-\frac{1+{i}}{\sqrt{2}}\zeta \right) \right\rbrace \right]\,, \end{equation} \tag{ 46 }$$

where ζ = exp ((z₀ − z)/δ) and I₁ is the integral

$$\begin{equation} I_1(\zeta) = \int _0^\infty \!{e}^{-\zeta s}\frac{s\,{d}s}{1+s^4}\,. \end{equation} \tag{ 47 }$$

The weakly magnetic and magnetically dominated regions correspond to ζ ≪ 1 and ζ ≫ 1, respectively.

Figure 5 shows the vertical profiles of $\hat{u}_z$ and $\hat{u}_x$, in the vicinity of the magnetic transition layer, from the same numerical solution shown in Figure 3. Also shown are their analytical counterparts, Equation (46), and the corresponding analytical profile of $\hat{u}_x$.

Zoom In Zoom Out Reset image size

Within the magnetically dominated region, where ζ ≫ 1, the flow is exponentially weak. From Equation (46), we find that

$$\begin{equation} \hat{u}_z\sim \frac{{i}k\delta u_{\rm t}}{\zeta ^{2}} = {i}k\delta u_{\rm t}\exp \left(\frac{2(z-z_0)}{\delta }\right). \end{equation} \tag{ 48 }$$

Within the weakly magnetic region, where ζ ≪ 1, $\hat{u}_z$ tends to a constant, w_t say, as expected from the Taylor–Proudman constraint. Specifically,

$$\begin{equation} w_{\rm t} = \hat{u}_z|_{\zeta =0} = -{i}\frac{\pi }{4}k\delta u_{\rm t}. \end{equation} \tag{ 49 }$$

The vertical mass flux within the weakly magnetic region is therefore tied to the differential rotation via the magnetic transition layer.

3.2.2. Global Solution

We now use relation (49) derived from our transition-layer solution to construct an approximate analytical solution for the global flow. In this way, we derive an explicit relation between the constant vertical velocity w_t within the weakly magnetic region and the prescribed differential rotation u_cz(z) within the convection zone. To construct the global solution, we must also match the flow within the weakly magnetic region, z₀ < z < h, to the flow within the convection zone, h < z < 1. The matching procedure is similar to that employed by GB10 and is described in Appendix B. Here we present only the result for w_t.

In the absence of stratification, n_rz = 0, we find that w_t is given by

$$\begin{equation} w_{\rm t} \;\; = \;\;\frac{{-{i}\bar{u}_{\rm cz}} }{\frac{kd}{2\tau _c} {\rm coth} \left(\frac{1-h}{d}\right) + \frac{4}{\pi k\delta } }, \end{equation} \tag{ 50 }$$

where d is the "convective" lengthscale⁴

$$\begin{equation} d = \frac{\sqrt{4\tau _c^2 + 1}}{k}, \end{equation} \tag{ 51 }$$

and where $\bar{u}_{\rm cz}$ is a weighted average of the forcing in the convection zone,

$$\begin{equation} \bar{u}_{\rm cz} \;\;=\;\; \frac{ \int _h^1\!\!{d}z\,u_{\rm cz}(z)\cosh \left(\frac{1-z}{d}\right)}{ \int _h^1\!\!{d}z\,\cosh \left(\frac{1-z}{d}\right)}\;. \end{equation} \tag{ 52 }$$

Since d > 1/k > (1 − h), the weight distribution in Equation (52) is roughly uniform, and so $\bar{u}_{\rm cz}$ is numerically close to the vertical average of u_cz(z) over the convection zone.

Figure 6 shows the vertical profiles of $\hat{u}_z$ in a series of numerical solutions with various magnetic scale heights δ. The bottom panel in Figure 6 compares the value of $\hat{u}_z$ in the weakly magnetic region to the value for w_t predicted by Equation (50). The fit is excellent provided that the magnetic transition layer is located sufficiently far above the lower boundary, z = 0, and sufficiently far below the base of the convection zone, z = h. Using Equation (44), these two conditions can be expressed as lower and upper bounds on Λ:

$$\begin{equation} \frac{2}{(k\delta)^2} \, {\ll} \, \Lambda \,{\ll} \, \frac{2{e}^{2h/\delta }}{(k\delta)^2}. \end{equation} \tag{ 53 }$$

Zoom In Zoom Out Reset image size

We have found that a confined magnetic field is able to halt the burrowing of the convection zone's differential rotation into the radiative interior, provided that the magnetic field strength exceeds a critical value B_crit given by Equation (45). Because the field strength in our model increases monotonically with depth, the radiative zone divides into an upper, weakly magnetic region, wherein B₀(z) ≪ B_crit, and a lower, magnetically dominated region, wherein B₀(z) ≫ B_crit. In the weakly magnetic region, the flow is subject to the Taylor–Proudman constraint, and the differential rotation u_x is independent of z. In the magnetically dominated region, the flow is subject to the Ferraro constraint, which suppresses both differential rotation and meridional flows. In our model B₀(z) increases exponentially with depth on a fixed lengthscale δ, and so the transition between the weakly magnetic and magnetically dominated regions occurs across a layer of thickness ≃ δ, within which B₀ ≃ B_crit.

The transition layer regulates the mass flux that downwells from the convection zone, as expressed by Equation (49). This is closely related to the mechanical constraint mentioned in Section 1 and can be explained as follows. The latitudinal differential rotation below the base of the convection zone extends downward until it meets the transition layer, where it winds up the magnetic field lines, producing a Lorentz torque that overcomes Taylor–Proudman balance and gyroscopically pumps a meridional flow. In the steady state, the downwelling mass flux within the weakly magnetic region must be exactly that demanded by the transition layer, which is expressed by Equation (49). An analogy can be made with the more familiar Ekman pumping that occurs in Ekman layers. Indeed, if the transition layer were replaced by an artificial frictional, impenetrable, and non-magnetic horizontal boundary, then Equation (49) would be replaced by the Ekman-pumping formula,

$$\begin{equation} w_{\rm t} = -\frac{1}{2}{i}k\delta _{\rm E}\, u_{\rm t} \end{equation} \tag{ 54 }$$

(GAA09), where δ_E = E^1/2_ν is the non-dimensional Ekman-layer thickness. In the solutions presented here we have δ ≫ δ_E, and so the magnetic transition layer pumps much stronger meridional flows than would be pumped by an artificial Ekman layer.

The solutions we have found bear many similarities to the tachocline model originally proposed by Gough & McIntyre (1998). In particular, we have a differentially rotating region below the base of the convection zone, which we might call the tachocline, and a thin magnetic boundary layer at the base of this region, which Gough & McIntyre called the "tachopause." The entire region below the tachopause is held in uniform rotation by the confined magnetic field. However, an obvious shortcoming of our unstratified solutions is that the "tachocline" has no vertical shear, quite unlike the solar tachocline. In the Gough & McIntyre model, the tachocline's vertical shear is explained by the presence of stable stratification via the thermal-wind relation. In the following sections, we reintroduce stratification into our model in order to better approximate conditions within the solar tachocline.

We begin again by exploring the behavior of the solutions of Equations (22)–(29) with varying governing parameters. In Figure 7, we show the vertical profiles of $\hat{u}_z$ in a series of numerical solutions with increasing stratification, measured in terms of σ. When σ is sufficiently small, ≲ 10⁻³, the profiles are indistinguishable from one another and from the unstratified solution, which has σ = 0. The same is true for the profiles of $\hat{u}_x$ in these solutions, which are shown in Figure 8(a). In particular, the weakly magnetic region (z₀ < z < h) remains in Taylor–Proudman balance. We refer to all such cases as "unstratified," since any effects due to stratification are negligible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As σ is increased beyond 10⁻³, we observe a reduction in $\hat{u}_z$ in the weakly magnetic region (see Figure 7). At the same time, the jump in $\hat{u}_x$ across the magnetic transition layer drops almost to zero, and instead $\hat{u}_x(z)$ rises smoothly and monotonically between z₀ and h, indicating that thermal-wind balance has been established within the weakly magnetic region. For brevity, and because of the obvious similarity with the model of Gough & McIntyre (1998), from here on we refer to the weakly magnetic region and magnetic transition layer in such cases as the "tachocline" and "tachopause," respectively.

We find that cases with significant stratification can be further categorized as being either "weakly stratified," "moderately stratified," or "strongly stratified." These three regimes are now described in turn.

4.1.1. The Strongly Stratified Regime

4.1.2. The Weakly Stratified and Moderately Stratified Regimes

Figure 8(b) shows vertical profiles of $\hat{T}$ from the same series of solutions shown in Figure 7. In all cases with σ ≲ 10⁻¹ we find that there is little variation in the temperature field $\hat{T}$ or its vertical gradient $\partial \hat{T}/\partial z$ across the tachopause. However, in cases with σ ≳ 10⁻¹ we observe a significant change in the temperature gradient at z = z₀, indicating that the flow within the tachopause contributes to the global-scale thermal equilibrium. We refer to cases with 10⁻³ ≲ σ ≲ 10⁻¹ as "weakly stratified" and cases with 10⁻¹ ≲ σ ≲ 1 as "moderately stratified."

Figure 9 illustrates how the mass flux in the tachocline varies with the strength of the interior magnetic field in the weakly and moderately stratified regimes. As in the unstratified case, we find that the position of the magnetic transition layer moves upward as Λ increases, in a manner which is more-or-less independent of σ. However, in contrast with the unstratified case, the strength of the meridional flow within the tachocline varies with Λ (cf. Figure 4). For increasing values of Λ, the meridional flow becomes stronger as the tachocline becomes thinner. This result has implications for the solar tachocline, as will be discussed in Section 4.2.2.

Zoom In Zoom Out Reset image size

Figure 10 shows meridional cross-sections through two solutions in the weakly and strongly stratified regimes. As mentioned in Section 1, the solar tachocline has σ ≃ 0.2–0.4 (GAA09) and therefore is not strongly stratified in the sense used here. We anticipate that either the weakly stratified or the moderately stratified regime offers the best approximation to the dynamics of the tachocline. In the following sections, we analyze these two regimes in detail in order to quantify more precisely the parameter ranges over which they apply.

Zoom In Zoom Out Reset image size

We proceed, as in Section 3.2, by deriving an analytical boundary-layer solution for the tachopause, which can then be used to construct a global analytical solution. In order to simplify the analysis, we will neglect the viscous terms in the equations, and so our analytical solutions cannot be applied in the strongly stratified regime.

4.2.1. The Tachopause

In the unstratified, weakly stratified and moderately stratified solutions presented in Figure 8(b), the temperature field $\hat{T}(z)$ is roughly constant in the region z − z₀ = O(δ). (Recall that the moderately stratified cases are characterized by a jump in ${d}\hat{T}/{d}z$, rather than $\hat{T}$.) We can therefore analyze the structure of the tachopause in all three regimes simultaneously by approximating $\hat{T}$ as a constant, T_t. As in Section 3.2.1, we neglect viscosity and make the boundary-layer approximation (34). We then recover Equations (36)–(42) exactly, except that Equation (38) is replaced by

$$\begin{equation} \frac{{d}\hat{p}}{{d}z} - T_{\rm t} \:=\: \frac{\Lambda }{\delta }\hat{b}_y{e}^{-z/\delta } - \Lambda \frac{{d}\hat{b}_y}{{d}z}{e}^{-z/\delta }\,. \end{equation} \tag{ 56 }$$

Within the tachocline, where z − z₀ ≫ δ, all the terms involving Λ are negligible, and so we have

$$\begin{eqnarray} \hat{u}_y&\simeq & 0 \end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray} 2\hat{u}_x&\simeq & {-}{i}k\hat{p} \end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray} \frac{{d}\hat{p}}{{d}z} - T_{\rm t} &\simeq & 0 \end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray} {i}k\hat{u}_y+ \frac{{d}\hat{u}_z}{{d}z} &\simeq & 0\,. \end{eqnarray} \tag{ 60 }$$

We deduce that

$$\begin{equation} \left. \begin{array}{@{}rcl@{}}\hat{u}_x&\sim & u_{\rm t} - \frac{1}{2}{i}kT_{\rm t}(z-z_0) \\ \hat{u}_y&\rightarrow & 0 \\ \hat{u}_z&\rightarrow & w_{\rm t} \end{array} \right\rbrace \;\;\; \mbox{as} \;\;\; \frac{(z-z_0)}{\delta } \rightarrow +\infty \,, \end{equation} \tag{ 61 }$$

where the constants u_t and w_t represent the values of $\hat{u}_x$ and $\hat{u}_z$ immediately above the tachopause, or equivalently, at the bottom of the tachocline.

As in Section 3.2.1, we can combine all of our boundary-layer equations into a single equation for $\hat{u}_z$:

$$\begin{equation} \left(\frac{{d}}{{d}z}-\frac{2}{\delta }\right) \left[ \left({e}^{z/\delta }\frac{{d}}{{d}z}\right)^4 + \frac{1}{4}k^4\Lambda ^2 \right]\hat{u}_z= \frac{k^4\Lambda }{4\delta }T_{\rm t}\,{e}^{2z/\delta }\,. \end{equation} \tag{ 62 }$$

In Appendix A, we show that the unique solution of Equation (62) that satisfies Equation (61) is

$$\begin{eqnarray} \hat{u}_z&=& {i}k\delta u_{\rm t}\left[\raisebox {0.5cm}{} I_1(\zeta) - \frac{\pi }{2}\,{\rm Re}\left\lbrace \exp \left(-\frac{1+{i}}{\sqrt{2}}\zeta \right) \right\rbrace \right] \nonumber \\ &+&\! \frac{1}{2}k^2\delta ^2 T_{\rm t}\left[\raisebox {0.5cm}{} I_2(\zeta) - \frac{\pi }{2}\,{\rm Re}\left\lbrace \left(\gamma -{i}\frac{\pi }{4}\right)\exp \left(-\frac{1+{i}}{\sqrt{2}}\zeta \right) \right\rbrace \right],\nonumber\\ \end{eqnarray} \tag{ 63 }$$

where γ is the Euler–Mascheroni constant, γ  =  0.577..., ζ = exp ((z₀ − z)/δ), as in Section 3.2.1, and

$$\begin{eqnarray} I_1(\zeta) &= \int _0^\infty \!\frac{s\,{d}s}{{e}^{\zeta s}}\,\frac{1}{1+s^4}\,, \end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray} I_2(\zeta) &= \int _0^\infty \!\frac{s\,{d}s}{{e}^{\zeta s}}\,\frac{\gamma +\ln s}{1+s^4}\,. \end{eqnarray} \tag{ 65 }$$

Figure 11 shows the vertical profiles of $\hat{u}_z$ and $\hat{u}_x$, and their analytical counterparts, from the same numerical solution shown in the top panel of Figure 10.

Zoom In Zoom Out Reset image size

From Equation (63), we deduce a relation between the constants T_t, u_t, and w_t,

$$\begin{equation} w_{\rm t} = \hat{u}_z|_{\zeta =0} = -\frac{{i}\pi }{4}k\delta u_{\rm t} - \frac{\gamma \pi }{8}k^2\delta ^2 T_{\rm t}\,. \end{equation} \tag{ 66 }$$

It can be verified that Equation (66) holds, to reasonable accuracy, in each of the solutions presented in Section 4.1. Furthermore, in each case the term in Equation (66) involving T_t is found to be negligible, and so Equation (66) reduces to relation (49) derived for the unstratified transition layer. The vertical mass flux within the tachocline is therefore tied to the differential rotation much as in the unstratified case described in Section 3.

We can also use Equation (63) to quantify the role of the tachopause in the global-scale heat flow. Within the tachopause the thermal energy equation (27) becomes

$$\begin{eqnarray} &&\frac{n_{\rm rz}^2}{{E}_\kappa }\hat{u}_z= \frac{{d}^2\hat{T}}{{d}z^2} \end{eqnarray} \tag{ 67 }$$

after making the boundary-layer approximation (34). Therefore, the change in the vertical temperature gradient across the tachopause is

$$\begin{eqnarray} \left[\frac{{d}\hat{T}}{{d}z}\right]_{z=z_0-\delta }^{z=z_0+\delta } \;\; &= \;\; \frac{n_{\rm rz}^2}{{E}_\kappa }\int _{z_0-\delta }^{z_0+\delta }\hat{u}_z\,{d}z \nonumber \\ \;\; &\sim \;\; \frac{n_{\rm rz}^2}{{E}_\kappa }w_{\rm t}\delta. \end{eqnarray} \tag{ 68 }$$

The transition between the weakly stratified and moderately stratified regimes occurs when Equation (68) is of comparable magnitude to the temperature gradient within the tachocline (see Figure 8(b)).

4.2.2. Global Solution

We can use relations (66) and (68) to construct an approximate analytical solution for the global flow. The details of the solution procedure are set out in Appendix B. We find that the vertical velocity w_t in the tachocline is

$$\begin{equation} w_{\rm t} = \frac{{-{i}\bar{u}_{\rm cz}}}{\frac{kd}{2\tau _c}{\rm coth} \left(\frac{1-h}{d}\right) + \frac{4}{\pi k\delta } + \frac{n_{\rm rz}^2}{2k^2{E}_\kappa }\,[G_1 + G_2 - G_3]}\,, \end{equation} \tag{ 69 }$$

which reduces to Equation (50) in the absence of stratification, n_rz = 0. The dimensionless "geometrical" factors G₁, G₂, and G₃ are given by Equations (B30)–(B32) in the weakly stratified regime and by Equations (B33)–(B35) in the moderately stratified regime. We refer to these factors as geometrical because of their dependence on the tachocline thickness, D = h − z₀. However, they also depend on f₀, the thermal diffusivity enhancement factor defined in Equation (20). Moreover, G₂ depends on the lengthscale d defined by Equation (51) and hence on the forcing timescale τ_c.

We can use Equation (69) to quantify the boundary between the unstratified and weakly stratified regimes identified in Section 4.1. This boundary is located where the term in the denominator of Equation (69) involving n_rz becomes as large as the term involving δ. In general, the geometrical factors are of order unity or smaller, and so the unstratified regime corresponds to

$$\begin{eqnarray} &&\frac{n_{\rm rz}^2}{{E}_\kappa } \lesssim k/\delta \,. \end{eqnarray} \tag{ 70 }$$

For the solutions in Figures 7 and 8, which have E_ν = 10⁻⁸, k = 2, and δ = 0.02, this condition is equivalent to σ ≲ 10⁻³, in good agreement with the numerical results.

The boundary between the weakly and moderately stratified regimes occurs where the change in the temperature gradient across the tachopause, given by Equation (68), becomes comparable to the temperature gradient within the tachocline. As described in Appendix B, this condition can be expressed approximately as

$$\begin{eqnarray} &&\frac{n_{\rm rz}^2}{{E}_\kappa } \simeq 1/(k\delta ^3)\,. \end{eqnarray} \tag{ 71 }$$

For the solutions in Figures 7 and 8, this condition is equivalent to σ ≃ 0.025, which is consistent with (although somewhat smaller than) the value σ ≃ 0.1 found in the numerical solutions.

The bottom panels of Figures 7 and 9 verify that the analytical prediction (69) for the vertical flow velocity in the tachocline agrees with the numerical results described earlier. Since viscosity was neglected in the derivation of Equation (69), this good agreement demonstrates that the weakly stratified and moderately stratified regimes are both inviscid, that is, viscous forces do not play a significant role in the dynamics. In particular, the solutions do not have an Ekman layer at the radiative–convective interface, z = h. We believe that the existence of such an Ekman layer in some previous studies (e.g., Gilman & Miesch 2004; Rüdiger et al. 2005; Garaud & Brummell 2008) arises from the treatment of the interface in those studies. In all the results presented here, the change in buoyancy frequency n at the interface is smoothed over a lengthscale Δ, representing the depth of convective overshoot, that greatly exceeds the Ekman length δ_E = E^1/2_ν. In the other studies just mentioned, the interface was modeled either as an upper boundary with an imposed differential rotation or by a change in n over a lengthscale ≪δ_E. Although the precise thickness of the Sun's overshoot layer is not known, it is surely thicker than the Ekman length (ν/Ω_☉)^1/2 ≃ 30 m, and so we argue that there is no Ekman layer in the solar tachocline.

It is useful to compare the expression for w_t given by Equation (69) to the corresponding expression derived by GB10 for meridional circulations pumped between the outer and inner convective zones of young lithium-dip stars (see their Equation (23)). As in their model, the magnitude of the circulation induced by the convective stresses is proportional to $\bar{u}_{\rm cz}$, a weighted average of the prescribed differential rotation in the convection zone. It is then moderated by whichever term in the denominator of Equation (69) is largest. The first term is always of order unity or larger, and so w_t can never exceed $\bar{u}_{\rm cz}$ in magnitude. The second term, involving δ, limits w_t to a flow rate which can be accommodated by the magnetic transition layer, as discussed in Section 3; this is an example of the "mechanical" constraint mentioned in Section 1. The third term, involving n_rz, represents the thermal constraint described in Section 1. This constraint limits w_t to a flow rate for which the temperature perturbations created by the advection of the background stratification can be balanced by thermal diffusion, thereby maintaining thermal equilibrium.

The expression for w_t given by Equation (69) can be simplified considerably in the parameter regime appropriate to the solar interior. In the solar tachocline we have n²_rz/E_κ ≃ 8 × 10¹³, and so the tachocline is not in the unstratified regime according to the criterion (70) unless the tachopause is unrealistically thin, δ ≲ 10⁻¹³. Therefore, the denominator of Equation (69) is dominated by the thermal term. The weakly stratified and moderately stratified regimes, according to Equation (71), correspond to δ ≲ 10⁻⁵ and δ ≳ 10⁻⁵, respectively. The solar tachocline could plausibly be in either of these regimes. Fortunately, it can be shown that G₁, G₂, and G₃ follow similar scaling laws in either case, differing only by factors of order unity. We can therefore consider both regimes simultaneously.

We expect the tachopause to be much thinner than the tachocline, δ ≪ D, which is itself much thinner than the solar radius, D ≪ 1. Under these conditions it can be shown that

$$\begin{equation} G_1/G_2 \sim f_0(kD)^2 \end{equation} \tag{ 72 }$$

$$\begin{equation} \mbox{and}\quad G_1/G_3 \sim D/\delta \gg 1 \end{equation} \tag{ 73 }$$

in both the weakly and moderately stratified regimes. The right-hand side of Equation (72) corresponds to the ratio of the thermal adjustment timescales in the tachocline and convection zone, respectively. In the solar tachocline, thermal adjustment by radiative diffusion has a timescale of about 10⁴ years, whereas the thermal adjustment timescale in the convection zone, estimated using mixing-length theory, is about 1 year. We therefore assume f₀(kD)² ≫ 1, in which case G₁ ≫ G₂, G₃. This same assumption was made implicitly in the models of Spiegel & Zahn (1992) and Gough & McIntyre (1998), in both of which the top of the tachocline was assumed to be isothermal. Under this assumption we find that G₁ ∼ (kD)³, and so Equation (69) becomes

$$\begin{equation} w_{\rm t} \,{\sim}\, -{i}\bar{u}_{\rm cz} \left(\frac{2{E}_\kappa }{n_{\rm rz}^2 kD^3}\right)\,. \end{equation} \tag{ 74 }$$

We note that w_t then depends only indirectly on the structure of the interior magnetic field via the tachocline thickness D. The dependence of w_t on D can be understood physically as a consequence of thermal equilibrium and thermal-wind balance, as first discussed by Gough & McIntyre (1998) (see also McIntyre 2007, p. 194). Since the vertical shear across the tachocline is of order $\bar{u}_{\rm cz}/D$, thermal-wind balance implies that there must be a temperature perturbation $\hat{T}$ of order $\bar{u}_{\rm cz}/(kD)$. Such a perturbation cannot persist within the convection zone, where heat is transported very efficiently by convective motions, but can persist within the tachocline, where heat transport is less efficient. In the tachocline, temperature perturbations diffuse at the rate E_κ/D², and thermal equilibrium therefore requires that $n^2_{\rm rz}w_{\rm t} \sim {E}_\kappa \hat{T}/D^2 \sim ({E}_\kappa /D^2)(\bar{u}_{\rm cz}/(kD))$, from which Equation (74) follows immediately. Using order-of-magnitude estimates for the tachocline thickness, D ≃ 0.01, and differential rotation, $\bar{u}_{\rm cz} \simeq 0.1$, as well as n²_rz/E_κ ≃ 8 × 10¹³, we find from Equation (74) that w_t  ∼  10⁻⁹, dimensionally w_t  ∼  10⁻⁴ cm s⁻¹. We note, however, that this result is rather sensitive to the tachocline thickness D, which is not well constrained by observations.

Our model therefore recovers Gough & McIntyre's scaling for the meridional flow within the tachocline and goes further by clarifying the physical assumptions under which that scaling is derived and by identifying the necessary conditions for those assumptions to hold. Under more general conditions, the meridional flow strength w_t is given by Equation (69), and Equation (74) should be regarded instead as an approximate upper bound on w_t. The more general formula (69) may be of relevance to the interiors of solar-type stars that are more weakly stratified or that have thicker tachoclines.

Of course, we must be careful when applying the results of our idealized model, with its artificially confined magnetic field, to the solar interior, or indeed to the interiors of other stars. However, since the formula for w_t given by Equation (74) has no explicit dependence on the structure of the background magnetic field B₀, we argue that this result should hold for more general, less artificial field configurations than that considered here. Moreover, the physical processes that give rise to Equation (74) will be present in any self-consistent model of the solar interior, and so Equation (74) ought to be a robust result under solar-like conditions.

We now ask whether the meridional flows predicted by Equation (74) would be of sufficient magnitude to confine an interior magnetic field against outward diffusion. In our idealized model, the interior field B₀ is confined by an artificially imposed downwelling U₀. A first approximation to the nonlinear magnetic confinement problem can be obtained by setting U₀ = |w_t|. Using Equation (7), we can then relate the tachopause thickness δ to the tachocline thickness D:

$$\begin{eqnarray} \frac{{E}_\eta }{\delta } &\sim & \bar{u}_{\rm cz} \frac{2{E}_\kappa }{n_{\rm rz}^2 kD^3} \end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray} \Rightarrow \;\; \delta &\sim & \frac{{E}_\eta }{{E}_\kappa } \frac{n_{\rm rz}^2 k}{2\bar{u}_{\rm cz} } D^3 \,. \end{eqnarray} \tag{ 76 }$$

This expression is equivalent to Equation (7) in Gough & McIntyre (1998).

We can now predict the strength of the magnetic field within the tachopause, B_t say, using Equation (45). Here, our model differs significantly from that of Gough & McIntyre. They assumed that thermal-wind balance would hold within the tachopause, as well as in the bulk of the tachocline, and as a result the thickness of the tachopause in their model, δ_GM say, was tied to the strength of the stable stratification. They found that

$$\begin{equation} \delta _{\rm GM} = \left(\frac{2{E}_\kappa }{k^4n_{\rm rz}^2\Lambda _{\rm t}} \right)^{1/6}, \end{equation} \tag{ 77 }$$

where Λ_t = Λ_loc(z₀) = B²_t/(4πηρ₀Ω_☉) is the Elsasser number within the tachopause. In our model, on the other hand, the Lorentz force within the tachopause overcomes thermal-wind balance, and we find that

$$\begin{equation} \delta = \left(\frac{2}{k^2\Lambda _{\rm t}}\right)^{1/2} \end{equation} \tag{ 78 }$$

(see Section 3.2.1). Hence, whereas Gough & McIntyre find a very strong dependence of B_t on D, namely B_t∝D⁻⁹ (see their Equation (8)), we find by combining Equations (76) and (78) that B_t∝D⁻³, and more precisely,

$$\begin{equation} \frac{B_{\rm t}^2}{4\pi \rho _0 R_\odot ^2 \Omega _\odot ^2} \sim \frac{ 8\bar{u}^2_{\rm cz} }{ n_{\rm rz}^4 k^4 D^6 } \frac{E^2_\kappa }{{E}_\eta } \,. \end{equation} \tag{ 79 }$$

Taking k = 2, n²_rz/E_κ ≃ 8 × 10¹³, and E_η ≃ 3 × 10⁻¹⁴, we find

$$\begin{eqnarray} \frac{B_{\rm t}^2}{4\pi \rho _0 R_\odot ^2 \Omega _\odot ^2} &\sim & 3\times 10^{-5} \left(\frac{\bar{u}_{\rm cz}}{0.1}\right)^2 \left(\frac{D}{0.01}\right)^{-6} \end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray} \Rightarrow \;\; B_{\rm t} &\sim & 10^3 \left(\frac{\bar{u}_{\rm cz}}{0.1}\right) \left(\frac{D}{0.01}\right)^{-3}{\rm G}\,. \end{eqnarray} \tag{ 81 }$$

This estimate is rather higher than the ∼1 G field predicted by Gough & McIntyre. We note that torsional Alfvénic oscillations of a ∼1000 G field could explain the 1.3 year oscillation detected in the tachocline's angular velocity (Gough 2000).

In recent years there have been several attempts to achieve magnetic field confinement in self-consistent, nonlinear, global numerical models of the solar interior. Various approaches have been taken, including axisymmetric steady-state calculations (Garaud & Garaud 2008), as well as two-dimensional and three-dimensional time-dependent simulations (Garaud & Rogers 2007; Rogers & MacGregor 2011; Strugarek et al. 2011). However, none of these numerical models has so far obtained solutions that are satisfactorily close to solar observations, and none has achieved magnetic confinement over an extended range of latitudes. These failures have led some to conclude that the magnetic confinement picture of Gough & McIntyre is unworkable. However, the analytical results presented here suggest that magnetic confinement by meridional flows can be achieved in the parameter regime appropriate to the solar tachocline. We can also use our results (1) to explain the lack of magnetic confinement in existing numerical models and (2) to guide parameter selection for future models.

As shown in Section 4.1, the amplitude of meridional flows within the radiative interior decreases monotonically as the stratification is increased. Moreover, in a time-dependent model, the timescale for the burrowing of the flows, given by Equation (3), increases with stratification. In the strongly stratified regime, with σ ≳ 1, the burrowing of meridional flows is even slower than viscous diffusion. But if σ  <  1, as is the case in the tachocline, then the meridional flows are approximately inviscid. When modeling the solar interior, it is common practice to impose a rotation rate and stratification profile that are close to solar, in order to ensure a realistic separation between the dynamical timescales. The viscous, thermal, and magnetic diffusivities (ν, κ, η) are then made as small as possible, subject to computational constraints. However, adopting this strategy does not guarantee that σ < 1. Indeed, if the ratio of buoyancy and rotational frequencies n_rz is chosen to match the true tachocline value, then the condition σ < 1 requires the Prandtl number ν/κ to be very small,

$$\begin{equation} \nu /\kappa < \frac{1}{n_{\rm rz}^2} \simeq 10^{-5}\,. \end{equation} \tag{ 82 }$$

This condition is satisfied in the solar interior, but is beyond the reach of current numerical simulations, which instead typically have a Prandtl number much closer to unity, in order to reduce the numerical stiffness of the equations. This places such simulations in the strongly stratified regime in which (1) meridional flow velocity decays exponentially with depth below the convection zone and (2) viscous stresses make a significant dynamical contribution, even if the Ekman number is small, i.e., even if E_ν ≪ 1.

The numerical difficulty is, in fact, rather easily avoided by imposing a weaker stratification, and thus allowing for a larger Prandtl number. For example, if n_rz = 10 then Equation (82) requires only that ν/κ < 10⁻², which is readily achievable numerically. We then expect the meridional flow velocity in the tachocline to scale as Equation (74), provided that the radiative–convective interface remains approximately isothermal (see discussion below Equation (72)).

A further numerical difficultly is the need to resolve both the thin tachocline and the thinner tachopause. This difficultly can be alleviated by allowing the tachocline to be somewhat thicker than is observed in the Sun.⁵ According to Equation (76), to obtain a tachopause of thickness δ = 0.01 and a tachocline of thickness D = 0.1, for example, we require

$$\begin{equation} \frac{{E}_\eta }{{E}_\kappa }\frac{n_{\rm rz}^2 k}{2\bar{u}_{\rm cz}} \sim 10\,. \end{equation} \tag{ 83 }$$

If we choose n_rz = 10, as suggested above, as well as $\bar{u}_{\rm cz} \simeq 0.1$ and k = 2, then this requires a diffusivity ratio η/κ = E_η/E_κ ∼ 10⁻². The strength of the interior magnetic field must then be chosen in accordance with Equation (79), which requires

$$\begin{equation} \frac{B_{\rm t}^2}{4\pi \rho _0 R_\odot ^2 \Omega _\odot ^2} \sim 50{E}_\kappa \,. \end{equation} \tag{ 84 }$$

Finally, we must ensure that the Ekman length E^1/2_ν is smaller than δ = 0.01.

All of the constraints just described can be satisfied by choosing E_κ = 10⁻³, E_η = 10⁻⁵, and E_ν = 10⁻⁶, as well as n_rz = 10 and B²_t/(4πρ₀) = 0.05R²_☉Ω_☉². These parameter values should be numerically achievable; in fact, the suggested values of E_κ, E_η, and E_ν are very similar to those used by Brun & Zahn (2006).

B.1.1. Solution in the Convection Zone, z ∈ [h, 1]

In the convection zone, the governing equations are well approximated by

$$\begin{eqnarray} {-}2\hat{u}_y&= - \frac{\hat{u}_x- u_{\rm cz}(z)}{\tau _c} \end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray} 2\hat{u}_x&= -{i}k\hat{p} - \frac{\hat{u}_y}{\tau _c} \end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray} 0 &= -\frac{{d}\hat{p}}{{d}z} + \hat{T} - \frac{\hat{u}_z}{\tau _c} \end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray} 0 &= {i}k\hat{u}_y+ \frac{{d}\hat{u}_z}{{d}z} \end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray} 0 &= \frac{{d}^2\hat{T}}{{d}z^2} - k^2 \hat{T}\,. \end{eqnarray} \tag{ B5 }$$

The general solution for the temperature perturbation can be written as

$$\begin{equation} \hat{T} = a \cosh k(z-1) + b \sinh k(z-1), \end{equation} \tag{ B6 }$$

where a and b are integration constants. From the boundary condition T = 0 at z = 1, we deduce immediately that a = 0. Combining the remaining equations yields

$$\begin{equation} \left(\frac{{d}^2}{{d}z^2} - \frac{1}{d^2}\right)\hat{u}_z= \left[-\hat{T} + \frac{2{i}}{k}\frac{{d}u_{\rm cz}}{{d}z}\right]\frac{\tau _c}{d^2}, \end{equation} \tag{ B7 }$$

where d is the lengthscale defined in Equation (51). We write the general solution for $\hat{u}_z$ as

$$\begin{equation} \hat{u}_z= - \frac{2{i}\tau _c}{kd^2}\int _z^1\!\!{d}z^{\prime }\,u_{\rm cz}(z^{\prime })\cosh \left(\frac{z^{\prime }-z}{d}\right) - \frac{b}{4\tau _c}\sinh k(z-1) + A\cosh \left(\frac{z-1}{d}\right) + B\sinh \left(\frac{z-1}{d}\right), \end{equation} \tag{ B8 }$$

where A and B are two additional integration constants. The boundary condition $\hat{w}=0$ at z = 1 implies that A = 0. (Since we have neglected the viscous terms in Equations (B1)–(B3) we cannot impose the stress-free boundary condition at z = 1. Including the viscous terms would lead to an Ekman-type boundary layer forming at z = 1, but the effect on the solution within the bulk of the convection zone would be of order E_ν ≪ 1.)

Finally, the pressure perturbation is found to be

$$\begin{eqnarray} \hat{p} &= \frac{2{i}}{k}u_{\rm cz} - \frac{d^2}{\tau _c}\frac{{d}\hat{u}_z}{{d}z} \end{eqnarray} \tag{ B9 }$$

$$\begin{eqnarray} &= - \frac{2{i}}{kd}\int _z^1\!\!{d}z^{\prime }\,u_{\rm cz}(z^{\prime })\sinh \left(\frac{z^{\prime }-z}{d}\right) + \frac{bkd^2}{4\tau _c^2}\cosh k(z-1) - \frac{Bd}{\tau _c}\cosh \left(\frac{z-1}{d}\right)\,. \end{eqnarray} \tag{ B10 }$$

B.1.2. Solution in the Tachocline, z ∈ [z₀, h]

Within the tachocline, we have

$$\begin{eqnarray} -2\hat{u}_y&= 0 \end{eqnarray} \tag{ B11 }$$

$$\begin{eqnarray} 2\hat{u}_x&= -{i}k\hat{p} \end{eqnarray} \tag{ B12 }$$

$$\begin{eqnarray} 0 &= -\frac{{d}\hat{p}}{{d}z} + \hat{T} \end{eqnarray} \tag{ B13 }$$

$$\begin{eqnarray} 0 &= {i}k\hat{u}_y+\frac{{d}\hat{u}_z}{{d}z} \end{eqnarray} \tag{ B14 }$$

$$\begin{eqnarray} \frac{n_{\rm rz}^2}{{E}_\kappa }\hat{u}_z&= \frac{{d}^2\hat{T}}{{d}z^2} - k^2 \hat{T}\,. \end{eqnarray} \tag{ B15 }$$

Equations (B11) and (B14) together imply that $\hat{u}_z$ is a constant, $\hat{u}_z= w_{\rm t}$ say. The remaining equations can then be integrated, and the general solution can be written as

$$\begin{eqnarray*} \hat{u}_z&=& w_{\rm t} \nonumber \\ \hat{T} &=& \left(T_{\rm t} + \frac{n_{\rm rz}^2}{k^2{E}_\kappa } w_{\rm t}\right)\cosh k(z-z_0) + K \sinh k(z-z_0) - \frac{n_{\rm rz}^2}{k^2{E}_\kappa } w_{\rm t} \nonumber \\ \hat{u}_x&=& u_{\rm t} - \frac{1}{2}{i}\left(T_{\rm t} + \frac{n_{\rm rz}^2}{k^2{E}_\kappa } w_{\rm t}\right)\sinh k(z-z_0) - \frac{1}{2}{i}K(\cosh k(z-z_0) - 1) + \frac{{i}n_{\rm rz}^2w_{\rm t}}{2k {E}_\kappa }(z-z_0) \nonumber \\ \hat{p} &=& \frac{2{i}u_{\rm t}}{k} + \frac{1}{k}\left(T_{\rm t} + \frac{n_{\rm rz}^2}{k^2{E}_\kappa } w_{\rm t}\right)\sinh k(z-z_0) + \frac{K}{k}(\cosh k(z-z_0) - 1) - \frac{n_{\rm rz}^2w_{\rm t}}{k^2 {E}_\kappa }(z-z_0), \nonumber \end{eqnarray*}$$

where u_t and T_t are the values of $\hat{u}_x$ and $\hat{T}$ at the bottom of the tachocline, and K is an additional integration constant.

B.1.3. Solution in the Magnetically Dominated Region, z < z₀

In this region, the differential rotation and meridional flow both vanish. The temperature perturbation $\hat{T}$ therefore satisfies

$$\begin{eqnarray} 0 &= \frac{{d}^2\hat{T}}{{d}z^2} - k^2\hat{T}\qquad \end{eqnarray} \tag{ B16 }$$

$$\begin{eqnarray} \Rightarrow \hat{T} &= \alpha \cosh kz + \beta \sinh kz\,. \end{eqnarray} \tag{ B17 }$$

Since $\hat{T}=0$ at the lower boundary z = 0, we must have α = 0.

The values of the seven integration constants b, B, w_t, u_t, T_t, K, and β are now determined by applying matching conditions across the interfaces z = z₀ and z = h. At the radiative–convective interface, z = h, we impose that $\hat{u}_z$, $\hat{p}$, $\hat{T}$, and $f(z)\frac{{d}\hat{T}}{{d}z}$ are all continuous.⁶ These are the same continuity conditions imposed by GB10, except that we use a more realistic thermal energy equation (9), and as a result we require the continuity of the heat flux, rather than the temperature gradient. Across the magnetic transition layer, at z = z₀, we impose continuity of $\hat{T}$ and relations (A28) and (A32) derived from our transition-layer solution.

The matching conditions lead to the following seven relations between the integration constants:

$$\begin{eqnarray} {-}\frac{2{i}\tau _c}{kd^2}\int _h^1\!\!{d}z^{\prime }\,u_{\rm cz}(z^{\prime })&\cosh \left(\frac{z^{\prime }-h}{d}\right) + \frac{b}{4\tau _c}\sinh k(1-h) - B\sinh \left(\frac{1-h}{d}\right) = w_{\rm t} \end{eqnarray} \tag{ B18 }$$

$$\begin{eqnarray} &&{-}\frac{2{i}}{kd}\int _h^1\!\!{d}z^{\prime }\,u_{\rm cz}(z^{\prime })\sinh \left(\frac{z^{\prime }-h}{d}\right) + \frac{bkd^2}{4\tau _c^2}\cosh k(1-h) - \frac{Bd}{\tau _c}\cosh \left(\frac{1-h}{d}\right) \nonumber \\ &&\qquad = \frac{2{i}u_{\rm t}}{k} + \frac{1}{k}\left(T_{\rm t} + \frac{n_{\rm rz}^2}{k^2{E}_\kappa } w_{\rm t}\right)\sinh kD+ \frac{K}{k}(\cosh kD - 1) - \frac{n_{\rm rz}^2w_{\rm t}}{k^2{E}_\kappa }D \end{eqnarray} \tag{ B19 }$$

$$\begin{eqnarray} {-}b\sinh k(1-h) &= \left(T_{\rm t} + \frac{n_{\rm rz}^2}{k^2{E}_\kappa } w_{\rm t}\right)\cosh kD + K\sinh kD - \frac{n_{\rm rz}^2w_{\rm t}}{k^2{E}_\kappa } \end{eqnarray} \tag{ B20 }$$

$$\begin{eqnarray} bk\cosh k(1-h) &= \frac{k}{f_0}\left(T_{\rm t} + \frac{n_{\rm rz}^2}{k^2{E}_\kappa } w_{\rm t}\right)\sinh kD + \frac{k}{f_0}K\cosh kD \end{eqnarray} \tag{ B21 }$$

$$\begin{eqnarray} T_{\rm t} &= \beta \sinh kz_0 \end{eqnarray} \tag{ B22 }$$

$$\begin{eqnarray} Kk &= \beta k\cosh kz_0 + \frac{n_{\rm rz}^2\delta }{{E}_\kappa }\left(\left(\frac{\pi }{4}\right)^3\frac{1}{2}k^2\delta ^2 T_{\rm t} - \gamma w_{\rm t} \right) \end{eqnarray} \tag{ B23 }$$

$$\begin{eqnarray} w_{\rm t} &= -{i}\frac{\pi }{4}k\delta u_{\rm t} - \frac{\pi }{8}\gamma k^2\delta ^2 T_{\rm t}, \end{eqnarray} \tag{ B24 }$$

where D = h − z₀ is the tachocline thickness.

We now seek an explicit expression for the value of the vertical flow w_t within the tachocline. We begin by eliminating B between Equations (B18) and (B19), which leads to

$$\begin{eqnarray} \frac{{i}kd}{\tau _c}w_{\rm t}&\coth \left(\frac{1-h}{d}\right) + \frac{{i}n_{\rm rz}^2w_{\rm t}}{k^2{E}_\kappa }kD - {i}\left(T_{\rm t} + \frac{n_{\rm rz}^2}{k^2{E}_\kappa } w_{\rm t}\right)\sinh kD - {i}K(\cosh kD - 1) \nonumber \\ &= 2\bar{u}_{\rm cz} - 2u_{\rm t} + \frac{{i}bkd}{4\tau _c^2}\left[ \frac{\sinh k(1-h)}{\tan\,h \big(\frac{1-h}{d}\big)} - kd\cosh k(1-h) \right], \end{eqnarray} \tag{ B25 }$$

where $\bar{u}_{\rm cz}$ is the weighted average of the forcing in the convection zone defined by Equation (52). The remaining six integration constants, including w_t, therefore depend on u_cz(z) only through its average $\bar{u}_{\rm cz}$. Hence the entire solution below the convection zone is determined by $\bar{u}_{\rm cz}$.

Next, we combine Equations (B22) and (B23) to eliminate β. The result can then be combined with Equations (B20) and (B21) to express b, T_t, and K in terms of w_t. The general result is rather complicated, so for simplicity we describe only two limiting cases, which correspond to the "weakly stratified" and "moderately stratified" regimes identified in Section 4.1.2.

B.2.1. The Weakly Stratified Regime, n²_rz/E_κ ≪ 1/(kδ³)

In this regime, we find that the term in Equation (B23) involving n_rz becomes negligible. Since this term arises from the change in the temperature gradient across the transition layer, this regime corresponds to the "weakly stratified" regime described in Section 4.1.2. After neglecting this term, we find that

$$\begin{eqnarray} b \; &= \; \frac{n_{\rm rz}^2w_{\rm t}}{k^2{E}_\kappa } \frac{1 - {\rm sech\,}k(z_0+D) \cosh kz_0}{f_0\cosh k(1-h) \tanh k(z_0+D) + \sinh k(1-h)}\qquad\qquad\quad \end{eqnarray} \tag{ B26 }$$

$$\begin{eqnarray} T_{\rm t} \; &= \; -\frac{n_{\rm rz}^2w_{\rm t}}{k^2{E}_\kappa } \frac{f_0\tanh kz_0[1 - {\rm sech\,}kD] + \tanh k(1-h) \tanh kz_0 \tanh kD}{f_0[\tanh kz_0 + \tanh kD] + \tanh k(1-h)[\tanh kz_0 \tanh kD + 1]} \end{eqnarray} \tag{ B27 }$$

$$\begin{eqnarray} K \; &= \; -\frac{n_{\rm rz}^2w_{\rm t}}{k^2{E}_\kappa } \frac{f_0[1 - {\rm sech\,}kD] + \tanh k(1-h) \tanh kD}{f_0[\tanh kz_0 + \tanh kD] + \tanh k(1-h)[\tanh kz_0 \tanh kD + 1]}. \end{eqnarray} \tag{ B28 }$$

B.2.2. The Moderately Stratified Regime, n²_rz/E_κ ≫ 1/(kδ³)

In this regime, the transition-layer temperature T_t turns out to be smaller than the value given by Equation (B27) by a factor E_κ/(n²_rzkδ³) ≪ 1. To good approximation, therefore, we can neglect T_t in each of the matching conditions (B18)–(B24), which is equivalent to setting z₀ = 0 in Equations (B26)–(B28).

After eliminating u_t between Equations (B25) and (B24), and applying the formulae for b, T_t, and K just derived, we arrive at an explicit formula for w_t,

$$\begin{equation} w_{\rm t} \;\; = \;\; \frac{-{i}\bar{u}_{\rm cz}}{\frac{kd}{2\tau _c}\coth \left(\frac{1-h}{d}\right) + \frac{4}{\pi k\delta } + \frac{n_{\rm rz}^2}{2k^2{E}_\kappa }[G_1 + G_2 - G_3]}, \end{equation} \tag{ B29 }$$

where

$$\begin{eqnarray} \quad G_1 &= kD - \frac{f_0[2 - 2{\rm sech\,}kD + \tanh kz_0 \tanh kD] + \tanh k(1-h) \tanh kD}{f_0[\tanh kz_0 + \tanh kD] + \tanh k(1-h)[\tanh kz_0 \tanh kD + 1]} \end{eqnarray} \tag{ B30 }$$

$$\begin{eqnarray} \ G_2 &= \frac{kd}{k^2d^2-1}\left[ kd - \frac{\tanh k(1-h)}{\tanh \big(\frac{1-h}{d}\big)} \right] \frac{1 - {\rm sech\,}kh \cosh kz_0}{f_0\tanh kh + \tanh k(1-h)}\qquad \end{eqnarray} \tag{ B31 }$$

$$\begin{eqnarray} G_3 &= \gamma k\delta \frac{f_0[1 - {\rm sech\,}kD]\tanh kz_0 + \tanh k(1-h) \tanh kD \tanh kz_0}{f_0[\tanh kz_0 + \tanh kD] + \tanh k(1-h)[\tanh kz_0 \tanh kD + 1]} \end{eqnarray} \tag{ B32 }$$

in the weakly stratified regime, and

$$\begin{eqnarray} G_1 &= kD - \frac{f_0[2 - 2{\rm sech\,}kD] + \tanh k(1-h) \tanh kD}{f_0\tanh kD + \tanh k(1-h)}\qquad\qquad \end{eqnarray} \tag{ B33 }$$

$$\begin{eqnarray} G_2 &= \frac{kd}{k^2d^2-1}\left[ kd - \frac{\tanh k(1-h)}{\tanh \big(\frac{1-h}{d}\big)} \right] \frac{1 - {\rm sech\,}kD}{f_0\tanh kD + \tanh k(1-h)} \end{eqnarray} \tag{ B34 }$$

$$\begin{eqnarray} G_3 &\simeq 0 \end{eqnarray} \tag{ B35 }$$

in the moderately stratified regime.