© ESO, 2012

1. Introduction

Turbulent thermal convection in the solar convection zone plays a key role for the maintenance of large scale flows (differential rotation, meridional flow) and solar magnetic activity. The angular momentum transport of convection maintains the global mean flows. Global flows relate to the generation of the global magnetic field, i.e. the solar dynamo. Differential rotation bends the pre-existing poloidal field and generates the strong toroidal field (the Ω effect), and mean meridional flow transports the magnetic flux equatorward at the base of the convection zone (Choudhuri et al. 1995; Dikpati & Charbonneau 1999). The internal structures of the solar differential rotation and the meridional flow are revealed by the helioseismology (see review by Thompson et al. 2003). Some mean field studies have reproduced these global flows (Kichatinov & Rüdiger 1993; Küker & Stix 2001; Rempel 2005; Hotta & Yokoyama 2011). These studies, however, used some kinds of parameterization of the turbulent convection, i.e. turbulent viscosity and turbulent angular momentum transport. Thus, a self-consistent thorough understanding of global structure requires the detailed investigation of the turbulent thermal convection. Turbulent convection is also important when describing the magnetic field itself. The strength of the next solar maximum in the predictions of the mean field model depends significantly on the turbulent diffusivity (Dikpati & Gilman 2006; Choudhuri et al. 2007; Yeates et al. 2008). In addition, the turbulent diffusion has an important influence on the parity of solar global field, the strength of polar field and so on (Hotta & Yokoyama 2010a,b).

There have been numerous LES numerical simulations of the solar and stellar convection (Gilman 1977; Gilman & Miller 1981; Glatzmaier 1984; Miesch et al. 2000, 2006; Brown et al. 2008) and their magnetic fields (Gilman & Miller 1981; Brun et al. 2004; Brown et al. 2010, 2011). In these studies, the anelastic approximation is adopted to avoid any of the difficulties caused by the high speed of sound. At the base of the convection zone, the speed of sound is about 200  km  s^-1. In contrast, the speed of convection is thought to be 50 m  s^-1 (Stix 2004). The time step must therefore be shorter owing to the CFL condition in an explicit fully compressible method, even when we are interested in phenomena related to convection. In the anelastic approximation, the equation of continuity is treated as (1)where ρ₀ is the stratified background density and v denotes the velocity. The anelastic approximation assumes that the speed of sound is essentially infinite, resulting in an instantaneous adjustment of pressure to flow changes. This is achieved by solving an elliptic equation for the pressure, which filters out the propagation of sound waves. As a result, the time step is only limited by the much lower flow velocity. However, owing to the existence of an elliptic the anelastic approximation has a weak point. When performing the numerical calculation using parallel computing, we must employ the frequent global communication. At the present time, the efficiency of scaling in parallel computing is saturated with about 2000–3000 CPUs in solar global simulations with the pseudo-spectral method (Miesch, priv. comm.). Higher resolution, however, is needed to understand the precise mechanism of the angular momentum and energy transport by the turbulent convection and the behavior of magnetic field especially in thin magnetic flux tube.

In this paper, we test the validity of a different approach to circumvent the severe numerical time step constraints in low Mach number flows. We use a method in which the speed of sound is reduced artificially by transforming the equation of continuity to (see, e.g., Rempel 2005) (2)where t denotes the time. Using this equation, the effective speed of sound becomes ξ times smaller, but the dispersion relationship for sound waves remains otherwise unchanged (wave speed decreases equally for all wavelengths). Since this technique does not change the hyperbolic character of the underlying equations, the numerical treatment can remain fully explicit, thus does not require global communication in parallel computing. We refer to the technique as the Reduced Speed of Sound Technique (RSST). This technique was used previously by Rempel (2005, 2006) in mean field models of solar differential rotation and non-kinematic dynamos, which essentially solve the full set of time-dependent axisymmetric MHD equations. Those solutions were however restricted to the relaxation toward a stationary state or very slowly varying problems on the timescale of the solar cycle. Here we apply this approach to thermal convection, where the intrinsic timescales are substantially shorter. To this end, we study two- and three-dimensional (3D) convection, in particular the latter will be non-stationary and turbulent.

The detailed setting of test calculation is given in Sect. 2. The results of our calculations are given in Sect. 3. We summarize our paper and give discussion of the RSST in Sect. 4.

2. Model

In this section, we describe the detailed settting of this study.

2.1. Equations

The two- or three-dimensional equation of continuity, equation of motion, equation of energy, equation of state, are solved in Cartesian coordinate (x,z) or (x,y,z), where x and y denote the horizontal direction and z denotes the vertical direction. The basic assumptions underlying this study are as:

• 1.

  the time-independent hydrostatic reference state;

• 2.

  that the perturbations caused by thermal convection are small, i.e. ρ₁ ≪ ρ₀ and p₁ ≪ p₀, here ρ₀ and p₀ denote the reference state values, whereas ρ₁ and p₁ are the fluctuations of density and pressure, respectively. Thus a linearized equation of state is used as Eq. (6);

• 3.

  that the profile of the reference entropy s₀(z) is a steady state solution of the thermal diffusion equation ∇·(Kρ₀T₀∇s₀) = 0 with constant K.

The formulations are almost identical to those of Fan et al. (2003), and given by where T₀(z), and s₀(z) denote the reference temperature and entropy, respectively, and e_z denotes the unit vector along the z-direction, γ is the ratio of specific heats, with the value for an ideal gas being γ = 5/3, and s₁ denotes the fluctuation of entropy from reference atmosphere. We note that the entropy is normalized by a value of the specific heat capacity at constant volume c_v. The quantity g is the gravitational acceleration, which is assumed to be constant. The quantity Π denotes the viscous stress tensor (7)and ν and K denote the viscosity and thermal diffusivity, respectively, and ν and K are assumed to be constant throughout the simulation domain. We assume for the reference atmosphere a weakly superadiabatically stratified polytrope where ρ_r, p_r, T_r, H_r, and δ_r denote the values of ρ₀, p₀, T₀, H₀ (the pressure scale height), and δ (the non-dimensional superadiabaticity) at the bottom boundary z = 0. Since |δ| ≪ 1, the value m is nearly equal to the adiabatic value, meaning that m = 1/(γ − 1). The strength of the diffusive parameter ν and K is expressed in terms of the non-dimensional parameters: the Reynolds number Re ≡ v_cH_r/ν and the Prandtl number Pr ≡ ν/K, where the velocity unit v_c ≡ (8δ_rgH_r)^1/2. We note that in this paper the unit of time is H_r/v_c. In all calculations, we set Pr = 1.

2.2. Boundary conditions and numerical method

We solve Eqs. (3)–(6) numerically. At the horizontal boundaries (x = 0,  L_x and y = 0,  L_y), periodic boundary conditions are adopted for all variables. At the top and the bottom boundaries, impenetrative and stress-free boundary conditions are adopted for the velocities and the entropy is fixed to At both the top and bottom boundaries (z = 0 and z = L_z), we set p₁ in the ghost cells such that the right hand side of the z-component of Eq. (4) is zero at the boundary (which is between ghost cells and domain cells), where the ghost cells are the cells beyond the physical boundary.

We adopt the fourth-order space-centered difference for each derivative. The first spatial derivatives of quantity q is given by (18)where i denotes the index of the grid position along a particular spatial direction. The numerical solution of the system is advanced in time with an explicit fourth-order Runge-Kutta scheme. The system of partial equations can be written as (19)for U_n + 1, which is the value at t_n + 1 = (n + 1)Δt is calculated in four steps The longest allowed time step, Δt_max, is determined by the CFL criterion. When both advection and diffusion terms are included in calculation, the time step reads (24)where (25)and c_tot is the total wave speed (26)where the effective speed of sound is expressed as (27)The time step determined by the diffusion term is (28)and c_ad and c_v are safety factors of order unity.

Using δ_r = 1 × 10^-4, the original speed of sound is about 35v_c at the bottom and 12v_c at the top boundary, respectively. In all calculations, Δx ~ 2.3 × 10^-2H_r. Thus, if we use c_ad = 1 and c_v = 1, Δt_ad = 6.4 × 10^-4  H_r/v_c and Δt_v = 1.5 × 10^-1  H_r/v_c. For all the ξ values considered in this paper, the time step remains restricted by the (reduced) speed of sound, thus the calculation is about ξ times more efficient with RSST.

3. Results

3.1. Two-dimensional study

We carried out four two-dimensional (2D) calculations with the RSST and one calculation without approximation (case 1–5). The values of the free parameters are given in Table 1. The superadiabaticity (δ) is 1 × 10^-6 at the bottom and about 2 × 10^-5 at the top boundary. It is almost the same superadiabaticity value as at the base of the solar convection zone. Figure  1 shows the time-development of entropy. In the beginning, non-linear time-dependent convection occurred (top panel), which changed to a steady state at later times (bottom panel). For the steady state, we compared the root mean squared (rms) velocity with different ξ. The rms velocity was defined as (29)Figure  2 shows our results where the value of ξ is from 1 to 80. The effective Mach number is defined as . If we have larger ξ value than 80, we cannot reach a stationary state, since there are some shocks generated by supersonic convection. The discussion of unsteady convection is given in the next session of 3D calculations. Even though the Mach number reaches 0.6 using ξ = 80 (panel a), the horizontal and vertical rms velocities are almost identical to those calculated with ξ = 1 (without RSST). The ratio of the rms velocities for each ξ and ξ = 1 are shown in Fig. 2e. The difference is always a few percent. This result is unsurprising since in the stationary state the equation of continuity becomes (30)whose solution must not depend on the value of ξ. We note that the cell size is to some degree affected by the aspect ratio of the domain. This does not affect our conclusions since this influence is the same for all considered values of ξ. To confirm the robustness of our conclusions, we repeated this experiment with a wider domain (26.16H_r × 2.18H_r instead of 8.72H_r × 2.18H_r) in which we find four steady convection cells with about 10% different rms velocities. In addition we find a dependence on ξ that is similar to that shown in Fig. 2, i.e. the solutions show differences of only a few percent as long as ξ < 80. In this section, we confirm that the RSST is valid for the 2D stationary convection when the effective Mach number is less than 1, i.e. ξ < 80 with δ_r = 1 × 10^-6. If we use δ_r = 1 × 10^-4 (result is not shown), the criterion becomes ξ < 8 for the 2D calculation.

Fig. 1 Time-development of entropy in a 2D calculation with parameters of case 1. (Top panel) t = 20. (Bottom panel) t = 400.

Fig. 2 Some quantities in 2D calculation. a) Maximum Mach number of each time step. b) Distribution of Mach number estimated with the rms velocity. c) Distribution of horizontal rms velocity vh. d) Distribution of vertical rms velocity. e) Ratio of the rms velocities for each ξ and ξ = 1.

As a next step, we investigate the dependence of the linear growth rate on ξ during the initial (time-dependent) relaxation phase toward the final stationary state. Figure  3 shows the linear growth of maximum perturbation density ρ₁ with different ξ. Black and red lines show the results with δ_r = 1 × 10^-6 and δ_r = 1 × 10^-4, respectively. These calculation parameters are not given in Table 1. In the calculation with δ = 1 × 10^-6, the growth rate decreases for values of ξ > 300, whereas it occurs at ξ > 30 in the calculation with δ_r = 1 × 10^-4. The reason can be explained as follows. In the convective instability, upflow (downflow) generates a positive (negative) entropy perturbation and then a negative (positive) density perturbation is generated by the sound wave. If the speed of sound is fairly slow, the generation mechanism of density perturbation is ineffective. In the calculation with δ_r = 1 × 10^-6 (δ_r = 1 × 10^-4), the growth rate with ξ = 200 (ξ = 20), however, is almost same as that with ξ = 1, even though flow with ξ = 200 (ξ = 20) is expected to be the Mach number Ma = 1.3, i.e. supersonic convection flow in the saturated state.

Fig. 3 Behavior in linear phase. Dependence of growth rate on ξ is shown. Black and red lines show the results with δr = 1 × 10-6 and 1 × 10-4, respectively.

3.2. Three-dimensional study

We now investigate the suitability of the RSST for 3D unsteady thermal convections (case 6–13). The value of superadiabaticity at the bottom boundary is 1 × 10^-4 and at the top 2 × 10^-3. Although this value is relatively large compared to solar value, the expected speed of convection is much smaller than the speed of sound, hence is small enough to allow us to investigate the validity of the RSST. The entropy of 3D convections with ξ = 1, 20, and 80 at t = 100H_r/v_c are shown in Fig. 4. The convection is completely unsteady and turbulent (animation is provided). The appearance of convection with ξ = 20 is almost the same as that with ξ = 1. We verify this using the Fourier transformation and auto-detection technique. However, the appearance of convection with ξ = 80 (bottom panel) differs completely from the others. This difference is most clearly visible in the animation of Fig. 4 that is provided with the online version.

Fig. 4 Snapshots of entropy of the 3D convection. Top, middle, and bottom panels correspond to ξ = 1, 10, and 80, respectively. Left (right) panels show entropy at top (bottom) boundary. (Animation is provided, and the difference between ξ = 1 and 80 is most clearly visible in the animation of Fig. 4, which is provided with the online version.)

Fig. 5 Some quantities averaged in time between t = 100 and 200Hr/vc in 3D calculations. a) Distribution of Mach number estimated with the rms velocity. b) Distribution of vertical rms velocity vh. c) Distribution of horizontal rms velocity. d) Distribution of rms power density of pressure and buoyancy. e) Distribution of rms power density of inertia.

We estimate the rms velocities for different ξ by calculating the average of values between t = 100 to 200H_r/v_c (Fig. 5: panel b and c). Without the RSST, i.e. ξ = 1 and using δ_r = 1 × 10^-4, the Mach number is 1 × 10^-2 at the bottom and 4 × 10^-2 at the top boundary. The rms velocities at ξ = 40 and 80 differ from those at the ξ = 1 by less than 15% and 30%, respectively. When we adopt ξ = 40 and 80, the Mach number estimated by rms velocity exceeds unity (Fig. 5: panel a), i.e. supersonic convection. This supersonic downflow frequently generates shocks and positive entropy perturbation, thus downflow is decelerated. This is the reason why the rms velocities with large ξ(=40, 80) are small. When ξ = 5, 10, 15, and 20, however, the rms velocities are in good agreement with that for ξ = 1. We estimate the rms power density of pressure, buoyancy, and inertia (Fig. 5). This results display almost the same trend as for the rms velocities. The power profiles ξ = 5, 10, 15, 20 agree with each other and those with ξ = 40, 80 disagree with those with ξ = 1. These results show that the RSST is a valid technique at least for ξ = 20, at which the Mach number is around 0.7. We note that the convection pattern among our 2D and 3D cases varies substantially. As a consequence, the ξ values for which the validity of RSST breaks down are also different in the 2D and 3D setups.

Fig. 6 Comparison of horizontal velocity spectra with different ξ. a) z = 2.18Hr b) z = 1.09Hr c) z = 0. Velocities are averaged in time between t = 100 and 200Hr/vc.

Fig. 7 Detection of convective cell. a) Original contour of density perturbation at z = 2.18Hr. b) Distribution of detected convection cell. Color and label (#n) show each convective cell.

In Fig. 6, we compare the average spectral amplitudes for different values of ξ. There is no significant difference between the spectral amplitudes for different values of ξ in the range from  1 to 20.

We investigate the distribution of cell sizes at the top boundary. This value is strongly related to the turbulent diffusivity and transport of either angular momentum or energy. The method for detecting the cell is as follows. At the boundary of the cell i.e., the region of downflow, the perturbation of density is positive and has a high value. When the density in a region exceeds a threshold, the region is regarded as the boundary of a convective cell. When a region is surrounded by one continuous boundary region, the region is defined as one convective cell. Figure  7 shows the detected cells. Each color and each label (#n) correspond to each detected cell. We estimate the size of all cells and compare the distribution of cell sizes for different ξ. The results are shown in Fig. 8. The cell size distribution follows a power law from about 0.01 to 10 and there is no dependence on ξ in the range from 1 to 20. Although when using this type of technique the size of cells tends to be large when neglecting smaller cells, our conclusion is not wrong, since all auto-detections are equally affected.

Fig. 8 Distribution of convective cell size with different ξ is shown. The cell size distribution is averaged in time between t = 100 and 200Hr/vc.

Fig. 9 The results with δr = 1 × 10-3. The format is the same as Fig. 5.

For the above two investigations, i.e. the Fourier analysis and study of cell detection, we conclude that the statistical features are unaffected by the RSST method as long as the effective Mach number (computed with the reduced speed of sound) does not exceed values of about 0.7, which corresponds to ξ = 20 in our setup.

To confirm our criterion that the RSST is valid if the Mach number is smaller than 0.7, we conduct calculations with larger superadiabaticity, i.e. δ_r = 1 × 10^-3 (cases 14–18). If our criterion is valid, required ξ with larger superadiabaticity must decrease. The results are shown in Fig. 9. Using δ_r = 1 × 10^-3, the calculations with ξ = 10, 15, and 20 generate a supersonic convection flow near the surface (Fig. 9a). It is clear that the results for ξ = 10, 15, and 20 differ from those with ξ = 1 and 5. The calculation with ξ = 5 corresponds to a flow whose Mach number is around 0.6. Thus, our criterion is not violated even for different values of the superadiabaticity.

Fig. 10 Dependence of  [∇·(ρ0v)] rms on ξ. In case 13, the value of ξ is 20 at the bottom and 4.5 at the top boundary.

Owing to our altered equation of continuity the primitive and conservative formulation of Eqs. (3) to (5) are no longer equivalent. For example, the equation of motion in conservative form is expressed as (31)where F denotes the pressure gradient, gravity, and Lorentz force. With some transformations, we can obtain (32)If the equation of continuity is satisfied, i.e. (∂ρ/∂t =  −∇·(ρv)), the primitive form is obtained as (33)However, using the RSST these two forms are no longer equivalent. We used the primitive formulation at the expense that energy and momentum are not strictly conserved; however, the consistency of our results for different values of ξ strongly indicates that this is not a serious problem for the setup we considered. Alternatively, we could also implement our modified equation of continuity in a conservative formulation. This would ensure that density, momentum, and energy are strictly conserved at the expense of a modified set of primitive equations. Figure 10 shows the dependence of  [∇·(ρ₀v)] _rms (=[ξ²∂ρ₁/∂t] _rms) on ξ. Using ξ = 5, this term is almost the same as that with ξ = 1. Although the deviation becomes large as ξ increases, it is not proportional to ξ².

In the previous discussion, we kept ξ constant in the entire computational domain. In the solar convection zone, the Mach number varies however substantially with depth, from  ~1 in the photosphere to  <10^-7 at the base of the convection zone, the speed of sound itself varies from about 7  km  s^-1 in the photosphere to 200  km  s^-1 at the base of the convection zone. A reduction in the speed of sound is therefore most significant in the deep convection zone, but not in the near surface layers. This could be achieved with a depth dependent ξ. Even if we use a conservative form of the equation of continuity, (34)the result must not be same as the result without any approximation of the statistical steady state. When the value averaged during a statistical steady state is expressed as  ⟨ a ⟩ , where a is a physical value, the equation of continuity becomes (35)with inhomogeneous ξ. The solution of Eq. (35) is different from the solution of the original equation of continuity in the statistical steady state, i.e. 0 = ∇·(ρ₀ ⟨ v ⟩ ). Thus, the statistical features such as rms velocity are not reproduced with an inhomogeneous ξ. If we use the non-conservative form of the equation of motion given in Eq. (2), this problem does not occur. Thus, we investigate the nonuniformity of ξ in case 13 using the non-conservative form, i.e., Eq. (3). To keep the Mach number uniform for all heights, we use ξ = 20/(δ(z)/δ_r)^1/2, i.e. ξ = 20 at the bottom and ξ = 4.5 at the top boundary. We note that the ratio of the speed of convection to the speed of sound is roughly estimated as . In this setting, the value of is 1 × 10^-2 at the bottom and 4 × 10^-2 at the top boundary. The same analysis as that for uniform ξ is done (see Figs. 5, 6, and 10). There is no significant difference between ξ = 20/(δ(z)/δ_r)^1/2 and ξ = 1. Although mass is not conserved locally for a non-homogeneous ξ using a non-conservative form of the equation of continuity, the horizontally averaged vertical mass flux is approximately zero in the statistically steady convection and the conservation total mass is not broken significantly. We conclude that an inhomogeneous ξ is valid under the previously obtained condition, i.e. the Mach number is less than 0.7.

4. Summary and discussion

We have applied RSST (described in Sect. 1) to 2D and 3D simulations of low Mach number thermal convection and confirmed the validity of this approach as long as the effective Mach number (computed with the reduced speed of sound) stays below 0.7 everywhere in the domain. The overall gain in computing efficiency that can be achieved therefore depends on the maximum Mach number that was present in the setup of the problem.

Since the Mach number is estimated to be 10^-4 in the base of the solar convection zone (Stix 2004), a several thousand times longer time step can be taken using the RSST. Therefore, the RSST and parallel computing with a large number of CPUs will make it possible to calculate large-scale solar convection with high resolution in the near future.

Compared to the anelastic approximation, there are three major advantages in RSST:

• 1.

  It can be easily implemented into any fully compressible code(regardless of the numerical scheme or grid structure) since itonly requires a minor change in the equation of continuity andleaves the hyperbolic structure of the equations unchanged.

• 2.

  Owing to the explicit treatment, it does not require any additional communication overhead, which makes it suitable for massive parallel computations.

• 3.

  The base of the convection zone and the surface of the sun where the anelastic approximation is broken can be connected, using a space-dependent ξ.

Overall we find that RSST is a very useful technique for studying low Mach number flows in stellar convection zones as it substantially alleviates any stringent time-step constraints without adding any computational overhead.

Acknowledgments

The authors want to thank Dr. M. Miesch for his helpful comments. Numerical computations were in part carried out on Cray XT4 at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. We would like to acknowledge high-performance computing support provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. The National Center for Astmospheric Reseach (NCAR) is sponsored by the National Science Foundation. This work was supported by the JSPS Institutional Program for Young Researcher Overseas Visits and the Research Fellowship from the JSPS for Young Scientists. We have greatly benefited from the proofreading/editing assistance from the GCOE program. We thank the referee for helpful suggestions that helped to improve our paper.

References

All Tables

All Figures

	Fig. 1 Time-development of entropy in a 2D calculation with parameters of case 1. (Top panel) t = 20. (Bottom panel) t = 400.
In the text

	Fig. 2 Some quantities in 2D calculation. a) Maximum Mach number of each time step. b) Distribution of Mach number estimated with the rms velocity. c) Distribution of horizontal rms velocity vh. d) Distribution of vertical rms velocity. e) Ratio of the rms velocities for each ξ and ξ = 1.
In the text

	Fig. 3 Behavior in linear phase. Dependence of growth rate on ξ is shown. Black and red lines show the results with δr = 1 × 10-6 and 1 × 10-4, respectively.
In the text

	Fig. 4 Snapshots of entropy of the 3D convection. Top, middle, and bottom panels correspond to ξ = 1, 10, and 80, respectively. Left (right) panels show entropy at top (bottom) boundary. (Animation is provided, and the difference between ξ = 1 and 80 is most clearly visible in the animation of Fig. 4, which is provided with the online version.)
In the text

	Fig. 5 Some quantities averaged in time between t = 100 and 200Hr/vc in 3D calculations. a) Distribution of Mach number estimated with the rms velocity. b) Distribution of vertical rms velocity vh. c) Distribution of horizontal rms velocity. d) Distribution of rms power density of pressure and buoyancy. e) Distribution of rms power density of inertia.
In the text

	Fig. 6 Comparison of horizontal velocity spectra with different ξ. a) z = 2.18Hr b) z = 1.09Hr c) z = 0. Velocities are averaged in time between t = 100 and 200Hr/vc.
In the text

	Fig. 7 Detection of convective cell. a) Original contour of density perturbation at z = 2.18Hr. b) Distribution of detected convection cell. Color and label (#n) show each convective cell.
In the text

	Fig. 8 Distribution of convective cell size with different ξ is shown. The cell size distribution is averaged in time between t = 100 and 200Hr/vc.
In the text

	Fig. 9 The results with δr = 1 × 10-3. The format is the same as Fig. 5.
In the text

	Fig. 10 Dependence of  [∇·(ρ0v)] rms on ξ. In case 13, the value of ξ is 20 at the bottom and 4.5 at the top boundary.
In the text