CC BY
16
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 18, ES1002, doi:10.2205/2018ES000615, 2018
Latitude dependence of convection and magnetic field generation in the cube
M. Yu. Reshetnyak1'2
Received 27 October 2017; accepted 16 January 2018; published 22 January 2018.
The 3D thermal convection in the Boussinesq approximation with heating from below and dynamo in the cube are considered. We study dependence of the convection intensity and magnetic field generation on the latitude in ,5-plane approximation. It is shown that kinetic energy gradually increases from the poles to the equator more than order of magnitude. The model predicts the strong azimuthal thermal wind, which direction depends on the sign of the thermal convective fluctuations. The spatial scale of the arising flow is comparable to the scale of the physical domain. The magnetic energy increases as well, however dynamo efficiency, i.e., the ratio of the magnetic energy to the kinetic one decreases to the equator. This effect can explain predominance of the dipole configuration of the magnetic field observed in the planets and stars. The approach is useful for modeling of the magnetohydrodynamic turbulence in planetary cores and stellar convective zones. KEYWORDS: Thermal convection; magnetic fields; planetary and stellar dynamo.
Citation: Reshetnyak, M. Yu. (2018), Latitude dependence of convection and magnetic field generation in the cube, Russ. J. Earth. Sci., 18, ES1002, doi:10.2205/2018ES000615.
1 Introduction energy sources in the model let to simulate the
various observable features of the magnetic fields, 3D modeling of convection and dynamo pro- like the spatial and temporal spectra of the mag-cesses in the planetary cores and convective zones netic fields, the butterfly diagrams in the case of of the stars is a modern branch of the magnetohy- the solar dynamo, reversals of the geomagnetic drodynamic (MHD) simulations [Riidiger at al., field.
2013]. There are two approaches which are usu- TT ...... , , ,. , .
However similarity of observations and simu-
ally used. For the large-scale fields modeling the
lations sometimes can be misleading because the
MHD equations are solved in the spherical geom- parameters used in 3D models are still quite far
etry, see, e.g.,[Jones, 2000]. Inclusion of the realis- from the desired ones. Briefly, the main problem
tic boundary conditions and distributions °f the is a turbulence which simulation requires resolu-
i^T ^TZi ^ r m, • r i.1, t i.u d tion of the small scales. The other difficulty is the xSchmidt Institute of Physics of the Earth of the Rus- . , J ...
sian Academy of Sciences, Moscow, Russia amsotropy of the flow, concerned with the rapid
2Pushkov Institute of Terrestrial Magnetism, Iono- daily rotation, as inthe case of the planetary cores.
sphere and Radio Wave Propagation of the Russian The spherical models which have dense distribu-
Academy of Sciences, Moscow, Russia tion of the mesh grid points near the poles, require the small time step (because of the Courant
Copyright 2018 by the Geophysical Center RAS. condition), and therefore do not suit to the tur-
http://rjes.wdcb.ru/doi/2018ES000615-res.html bulence modeling. Meanwhile, simulations in the
ES1002
1 of 7
Cartesian geometry with the homogeneous grids could be very helpful. These simulations reproduce the cascade processes between the scales, are helpful for estimates of the turbulent coefficients, and demonstrate various remarkable properties of the MHD turbulence [Brandenburg and Subra-manian, 2005].
In spite of the fact that these two approaches were used for years, the flat layer dynamo simulations did not take into account the latitude dependence to the moment. In other words the angle between the angular rotation axis O and gravity g was set to zero [Jones and Roberts, 2000], [Buffett, 2003], [Cattaneo et al., 2003]. This choice of parameters corresponds to the geographic poles. Having in mind that the most common configuration of the magnetic field in planets and stars is the dipole, and 3D spherical dynamo models predict existence of the toroidal counterpart, concentrated at the mean latitudes, it looks tempting to include the latitude dependence in the flat layer models immediately. The lack of the papers, devoted to this problem in the geodynamo and stellar applications, looks quite surprising because in the meteorology the latitude dependence, known as the ft -plane approximation, was already used in the Cartesian models for a long time [Pedlosky, 2012].
The other motivation of this paper is to distinguish the physical effects related to the angle between O and g definitely, leaving apart influence of the inner core and spherical boundaries, which can change hydrodynamics2, and the magnetic field generation substantially.
Below we consider the standard 3D MHD model in the rapidly rotating cube. The model includes the thermal convection equations in the Boussi-nesq approximation with the heating from below. The liquid is conductive and at the large magnetic Reynolds numbers the fluid motions can generate the magnetic field. To solve these equations we use MPI C++ code, based on the predictor-corrector method [Canuto et al., 1988]. We check how these processes depend on the angle between O and g. Some analytical predictions, based on the analogy with the motion of the charged particle in the electromagnetic field are considered as well.
2The sign of the boundaries curvature defines the direction of the Rossby waves propagation [Busse, 2002].
2. Dynamo in the Cube and Numerical Methods
The dimensionless dynamo equations for an incompressible fluid (V ■ V = 0) in the cube of the height L = 2n, rotating with the angular velocity O, in the Cartesian system of coordinates (x, y, z) have the form:
dA
— = V x B + q-1AA, B = rot A
dt
EPr
-1
dV
■w + (V v) V
rot B x B - V P — x V + RaT 1Z + EAV dT
— + (V -V) (T + To) = AT.
(1)
Velocity V, magnetic field B (derived from the vector potential A), pressure P and the typical diffusion time t are measured in units of k/L, a/2Q khp;, pK2/L2 and L2/ k respectively, where k is the thermal diffusivity, p is the density, ^
the permeability, Pr = — is the Prandtl number,
E = oot 2 is the Ekman number, v is the kine-2sZL
matic viscosity, ] is the magnetic diffusivity, and
q = k/t] is the Roberts number. Ra = L
2&2 k
is the modified Rayleigh number, a is the coefficient of the volume expansion, 6T is the unit of the temperature fluctuations T, is the gravitational acceleration, and T0 = 2tt — z is the temperature profile, corresponding to the heating from below.
The unit vector 1n defines direction of the angular velocity O. The angle between O and gravity g, directed along the z-axis, is equal to the colatitude , which is related to the latitude as d = 90° — 0.
The problem is closed with the periodical boundary conditions in the (x, y)-plane. In z-direction the following simplified boundary conditions
T = Vz = Az =
dVx dVy dAT, dA,
■ ■
■
■
0
at z = 0, 2^ are used. Conditions for A are the so-called pseudo-vacuum boundary conditions, correspond to the following conditions for the mag-
netic field: BT = R, =
dBz
0. Then the normal
y dz
component to the boundary of the electric current J = V x B is zero.
The system (1) was solved using the finite differences of the second-order in space and time. For approximation of the time derivative the three-layer time scheme was used:
df 3/n+1 - 4fn + fn-1
dt
2 St
3fn+1 - 4fn + f
•n I fn— 1
2 St
where
Fn = Cn + 1 V2/n,
3 V* - 4 Vn + V
n I \T'n—1
2 St
2 Fn — Fn-1 + 1V2V*,
(4)
where
Fn = - E Pr-1 (Vn -V) Vn + rot Bn x Bn-1n x Vn + RaTn 1Z + 2 E AVn.
(5)
Eqs(4),(5) lead to a Poisson-type equation for V*.
Then pressure Pn+1 was derived from the continuity equation by solving the another Poisson problem:
3
v2Pn+1 = — V • V* V 2 St V
with the Neumann boundary condition for z-coor-dinate:
dPn+1 3 _
-=-V *
dz 2 5t z
The last step provides incompressibility of the velocity field:
V
n+1
V* — — Vp n+1. 3
The second order up-wind scheme was used for approximations of the convective terms in the heat and the Navier-Stokes equations:
where n denotes the time step St.
For T and A the corresponding equations from (1) were written in the form:
V-—T- =
OX
= 2Fn — Fn-1 + 1V2 fn+1, 2 (2)
(3)
and С denotes the corresponding convective term. Eqs(2),(3) with respect to /ra+1 were solved using the Gauss-Seidel method.
While using the vector potential A provides divergence free of the magnetic field B, incompress-ibility of the velocity field V should be provided by some special technique. Here we use the predictor-corrector method [Canute et at., 1988], [Bandaru et at., 2016], introducing the intermediate velocity field V* by equation:
(3Tt - 4Ti-1 + Tl-2)Vl/(2Sx), V > 0 (-T+2 + 4Tm - 3Ti)Vi/(25x), V < 0,
(6)
where denotes the index of the grid step x. The scheme (6) was used for y-,z-directions in the same way.
This approach was realised in C++ code with MPI. The whole domain was divided into ( N x N), subdomains in (x, y) coordinates, where the MHD equations (1) were solved. The subdomains exchanged by its boundaries at the each time step n. In this paper the mesh grid (295 x 295 x 125) in (x, y, z) coordinates, and N = 6 were used. The simulations were done at two linux workstations Intel(R) Xeon(R) CPU E5-2640 with 40 cores and 128Hb common memory at the station. Each run required 2-4 days, depending on the time step , which was in the range (2 + 5) 10-7.
3. Pure Convection
Before to consider the full dynamo the pure turbulent convective regime with E = 4 10-5, Ra = 9103, Pr = 1 was studied. Ten runs with step 10° in the latitude § were performed. After some intermediate stage solutions reached the quasi-stationary states. To measure the intensity of convection we estimated the mean over the volume kinetic energy, see Figure 1, defined as
Ek=Ы v2(1 r3,
V
V = 8^3.
The considered regimes correspond to the developed turbulence with the Reynolds number Re = 2ityJ 2Ek/3 ^ 1. All the energies increase from the poles to the equator. Only E^ at the equator has sharp minimum. This behaviour is expectable because at the equator, § = 0, x-compo-nents of the Coriolis and Archimedean forces are
The latitude dependence of the kinetic denote the kinetic ener-
El, El
Figure 1
energies. E%, gies of Vx-, Vy-, and Vz-components of the velocity, respectively. The total kinetic energy E^ =
El + Ey + Ez,.
zero, and the non-zero value of Vx is provided by the non-linear term in the Navier-Stokes equation only. The ratio of the maximum and minimum of E^ is 22. This effect is quite strong and should be explained in some way.
Moreover, analysis of the spatial structure of the flows, see Figure 2, reveals that the small-scale cyclonic convection, existed at the poles, changed to the large-scale convection at the equator. Note, that the Reynolds number in the latter case is larger, and the flow is "more" turbulent, but in the same time it is large-scale. It should be noted that Vy-component at the equator is perpendicular to ft and it should be twisted at the small scale in the similar way, as it was at the poles.
To find the origin of the large-scale convection the analogy with the motion of the charged particle in the constant electromagnetic field, considered in the next section, is instructive.
4. Analogy with the Charged Particle Moving in the Electromagnetic Field
Similarly to the motion of the charged particle in the constant in time and homogeneous in space electromagnetic field, see [Arzimovich and Lukianov, 1972], we consider two limiting cases, where Archimedean force and angular rotation velocity of the system are directed along the same
axis, and the other, where these forces are perpendicular.
In the first, the most studied case, which corresponds to the geographic poles, the flow is accelerated by the Archimedean force ~ Ra T. Due to the Coriolice force any motion in the orthogonal plane to the gravity and ft is twisted with radius r, defined by relation Ro V± /r ~ V±, i.e. r ~ Ro V±, where V± is the velocity orthogonal to ft. The Rossby number Ro = Pr-1 E for the Earth's core ~ 10-15. Even with V± ~ Re = 109, estimated from the large-scale velocity, based on the west drift of the geomagnetic field, one has r ~ 10-6 in units of the liquid core's scale. Taking into account decrease of the kinetic energy spectrum will only decrease estimate of r. This estimate is valid for the large velocities with negligible viscous dissipation.
The linear analysis at the threshold of convection generation, where viscous diffusion is important, also predicts existence of the small scale in the perpendicular plane: r ~ E1/3 = 10-5 [Roberts, 1968], [Busse, 1970].
The both estimates demonstrate that the small-scale convection at the poles is a quite natural phenomenon, and it appears in the rapidly rotating objects with Ro = E • Pr ^ 1 even at the critical Rayleigh numbers.
For the other case, which corresponds to the equator plane, let Archimedean force is still directed along the z-axis and ft along x, and the initial velocity is at the (y, z)-plane. Then the trajectory of the particle remains in the same plane, and its motion is described by equations:
Ro y = z
Ro z = —y + Ra T,
(7)
where dot is for the time derivative.
Eqs (7), written in the reference system moving with the velocity v along the y-axis, after substitution y1 = y — vt, have the form:
Ro y\ = —z
Ro z = y1 + v + Ra T.
(8)
Choosing v = — Ra T, one has circular motion in (yi, z)-plane with frequency ~ Ro-1. The trajectory in the original (y, z)-plane is a trochoid, i.e. the superposition of the circular motion and the drift with velocity v. In terms of the spherical geometry y corresponds to the azimuthal direction. Addition of the initial velocity in x-direction,
Figure 2. The x-sections of ^-components (upper line) and z-sections of Vz-components (bottom) of the velocity field. The left column corresponds to the pole, and the right one - to the equator.
which remains constant because of absence of forces in this direction, does not change situation.
In (8) T is a fluctuation of the temperature relative to the non-convective distribution, and it can change the sign. By analogy with the motion of the charge in the magnetic field one has thermal separator, which divides the hot (T > 0, v < 0) and cold (T < 0,v > 0) flows.
Note that estimates of the radius of rotation in planes perpendicular to the axis of rotation coincide in the both cases. The difference is existence of a thermal wind with a velocity in the az-imuthal direction in the latter case. It is this wind has large spatial scale, already detected in Figure 2 at the equator. Due to the continuity equation the other components of the velocity can also posses the large-scale counterpart, as we observe it in Vz-flow.
5. Dynamo
Starting from the pure convection velocity and temperature distributions, obtained in Section 3, and the small magnetic field initial seed, the full dynamo system (1) were integrated in time up to the state where all the physical fields stabilised at the quasi-stationary regime. The corresponding estimates of the magnetic energy, defined as
E =
1
2Ro V
B2 d r3,
are ploted versus lattitude § in Figure 3.
The behaviour of the magnetic energy is similar to the kinetic energy up to some details near the equator plane. The largest increase of the magnetic energy with the latitude demonstrates
Figure 3. The latitude dependence of the magnetic energies.
By-component, which is stretched by the strong large-scale thermal wind Vy. In contrast to the total kinetic energy, Ek, the magnetic energy Em slightly decreases at the equator. In some sense it resembles behaviour of the magnetic field in the spherical models with the odd configurations of the magnetic field with respect to the equator, e.g., dipole. Increase of the magnetic energy Em relative to the poles can reach factor 4.
One can expect that the large-scale flow, based on the thermal wind, and the small-scale cyclonic convection have different efficiencies of the magnetic field generation. To test this hypothesis the
E
ratio £ = —m was plotted as a function of $ in Fig-
Ek
ure 4. It appears that in the range of 50° < $ < 90° £ is approximately constant and then decreases to the equator in one order of the magnitude. We conclude that the large-scale flow with the small gradients is less effective than the cyclonic convection with the non-zero net helicity [Reshetnyak, 2017].
6. Conclusions
Our simulations clearly demonstrate that the angle between the axis of rotation and gravity changes not only the magnitude of the mean parameters in the flat model, like energies, but the structure of the flow, its spectra, as well. These effects are already notable at the mean latitudes, where the magnetic energy maximum is localised in the spherical models, and should be taken into
account in future. Of course, due to variety of the physical effects, application of these results to the spherical shells, should be done carefully. Thus existence of the well-known in geodynamo since [Glatzmaier and Roberts, 1995], see also [Reshetnyak and Pavlov, 2016], the large-scaled vortexes in the Taylor cylinder at the large Rayleigh numbers, introduce the new complexity in the model.
The other issue is improvement of the numerical methods for the MHD turbulence modling. The finite difference methods are the most promising approach for the multi-core simulations. The modern higher order approximations of MHD differential equations are already comparable by accuracy to the spectral methods [Brandenburg and Subramanian, 2005], which for years were the best choice for such problems. The main advantage of the finite differences is their scalability at the supercomputers. We hope that this paper can be useful for the further development of the realistic MHD turbulence models with rotation. In spite of the incompressible form of convection's equations, considered above, the substitute of variables V ^ pV can be used for transition to the anelastic approximation, suitable to the dynamo in the compressible medium.
Acknowledgment. The author acknowledges financial support from Russian Science Foundation, the grant No16-17-10097.
Figure 4. The latitude dependence of the magnetic field generation efficiency
References
Arzimovich, L. A., S. Yu. Lukianov (1972), Motion of the charged particles in the electromagnetic fields, n/a pp., Nauka, M.. (in Russian) https://doi.org/10.1007 /978-3-642-84108-8 Bandaru, V., T. Boeck, D. Krasnov, J. Schumacher (2016), A hybrid finite difference/boundary element procedure for the simulation of turbulent MHD duct flow at finite magnetic Reynolds number, J. Comp. Phys., 304, No. 6, 320-339, https://doi.org/10.1016 /j.jcp.2015.10.007 Brandenburg, A., K. Subramanian (2005), As-trophysical magnetic fields and nonlinear dynamo, Phys. Rep., 417, No. 6, 1-209, https://doi.org/10. 1016/j.physrep.2005.06.005 Buffett, B. (2003), A comparison of subgrid-scale models for large-eddy simulations of convection in the Earth's core, Geophys. J. Int., 153, No. 3, 753-765, https://doi.org/10.1046/j.1365-246x.2003.01930.x Busse, F.-H. (1970), Thermal instabilities in rapidly rotating systems, Fluid Mech., 44, 441-460, https://doi.org/10.1017/S0022112070001921 Busse, F.H. (2002), Convective flows in rapidly rotating spheres and their dynamo action, Phys. Fluids, 14, 1301-1314, https://doi.org/10.1063/L1455626 Canuto, C., M. Y. Hussini, Q. A. Zang (1988), Spectral Methods in Fluids Dynamics, 31-75 pp., SpringerVerlag, Berlin. https://doi.org/10.1007/978-3-642-84108-8
Cattaneo, F., et al. (2003), On the interaction between convection and magnetic fields, Astrophys. J., 588, No. 2, 1183-1198, https://doi.org/10.1086/ 374313
Glatzmaier, G., et al. (1995), A three-dimensional convective dynamo solution with rotating and finitely
conducting inner core and mantle, Phys. Earth Planet. Int., 91, 63-75, https://doi.org/10.1016/0031-9201(95)03049-3 Jones, C. A. (2000), Convection-driven geodynamo models, Phil. Trans. R. Soc. London, A358, 873-897, https://doi.org/10.1098/rsta.2000.0565 Jones, C. A., P. H. Roberts (2000), Convection-driven dynamos in a rotating plane layer, J. Fluid Mech., 404, 311-343, https://doi.org/10.1017/S0022112099007 363
Pedlosky, J. (2012), Geophysical fluid dynamics,
999 pp., Springer-Verlag, Berlin. https://doi.org/10. 1007/978-1-4684-0071-7 Reshetnyak, M., V. Pavlov (2016), Evolution of the Dipole Geomagnetic Field. Observations and Models, Geomagnetism and Aeronomy, 56, No. 1, 110-124, https://doi.org/10.1134/S0016793215060122 Reshetnyak, M. (2017), The anisotropy of hydrody-namical and current helicity, Astronomy Reports, 61, No. 9, 783-790, https://doi.org/10.1134/S1063772 917080108
Roberts, P.-H. (1968), On the thermal instability of a rotating-fluid sphere containing heat sources, Phil. Trans. R. Soc, A263, No. 1, 93-117, https://doi.org/10.1098/rsta.1968.0007 Rudiger, G., R. Hollerbach, L. L. Kitchatinov (2013), Magnetic Processes in Astrophysics: Theory, Simulations, Experiments, 356 pp., Wiley-VCH, Weinheim, Germany. https://doi.org/10.1002/9783527648924
M. Yu. Reshetnyak, Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences, Moscow, Russia, Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation of the Russian Academy of Sciences, Moscow, Russia, ([email protected])
