non-dipolar, and failed dynamo cases, respectively. We categorized strong dipo-lar when the magnetic energy is dipo-larger than the kinetic energy in each simulation result, and vice versa. Sustaining the dynamo with a smaller inner core size re-quires a larger Rayleigh number. This is consistent with the findings of Heimpel et al. (2005). Atri/ro = 0.35, almost all the sustained dynamo cases were strong dipoles. Atri/ro = 0.25, there were strong dipolar dynamo cases and weak dipo-lar dynamo cases. At ri/ro = 0.15, there were weak dipolar and non-dipolar dynamo cases.
ri/ro = 0.25than the present Earth. In contrast, the dipolarity atri/ro = 0.15is smaller than the present Earth’s dipolarity in all cases. The dipole component is not dominant.
In numerical dynamos atri/ro = 0.35, we verified that the transition between the dipole and non-dipole isfdip ≈0.35(Christensen and Aubert, 2006; Olson et al., 2011). Our results are consistent with this transition. While dipolarity is an ef-fective index if dynamos can be categorized into large and small dipolarity groups, the combination of dipolarity and the ratio of extrapolation from fitting assesses the dipolar dominance if the dipolarity changes gradually, as in our results.
At ri/ro = 0.15, an axial dipole field formed by a single column (Heimpel et al., 2005). In this study, a dipole field is also formed by some azimuthally localized narrow columns around the dynamo-onset cases. Emag is found to be always smaller thanEkinin allRacases. The magnitude relationship is the same as that of Heimpel et al. (2005). A strong dipole is sustained with a smaller inner core in the fixed flux calculation (Hori et al., 2010), changing the core power based on the thermal history (Driscoll, 2016), or the buoyancy gained by light elements (Lhuillier et al., 2019). Clarifying how heat flow at boundaries sustains the dipole requires further numerical simulations.
Our proposed method of evaluating the dipolar dominance, fmag fit, enables quantitative investigations of the magnetic field structure in the past environment.
Knowledge from paleomagnetic analyses, such as VDM and VGP (virtual
ge-omagnetic pole), is acquired based on the assumption that the gege-omagnetic field was dipolar-dominated in the past (Merrill et al., 1996). In contrast, the VGP paths and actual behavior of the geomagnetic field are not dipolar-dominant. Investiga-tion of the numerical dynamo with our proposed method is capable of improving the understanding of the actual behavior of the geomagnetic field and paleomag-netic observations.
Table 3.1:Results ofEkin for thermal convection att = 6τν inri/ro=0.15,0.25, and 0.35under the FT boundary condition.
ri/ro Ra[×105] Ekin 0.15 1.0 8.31×10−6
0.15 1.2 5.72
0.15 1.25 8.15
0.15 1.3 10.62
0.15 1.35 13.16
0.15 1.4 15.80
0.15 1.45 18.54
0.25 0.70 8.55×10−4
0.25 0.75 2.41
0.25 0.78 4.94
0.25 0.80 6.62
0.25 0.82 8.35
0.25 0.85 11.02
0.25 0.90 15.70
0.35 0.55 1.76×10−4
0.35 0.58 2.38
0.35 0.60 4.54
0.35 0.62 6.76
0.35 0.65 10.22
0.35 0.67 12.60
0.35 0.70 16.28
Table 3.2:Results ofEkininri/ro=0.15,0.25, and0.35.
ri/ro Ra[×105] Ra/Racrit Ekin
0.15 760 7.0 758.7
0.15 870 8.0 926.7
0.15 980 9.0 1190
0.15 1100 10.1 1431
0.15 1300 11.9 1747
0.15 1500 13.8 2328
0.15 1700 15.6 2389
0.25 140 1.9 76.69
0.25 160 2.2 69.64
0.25 180 2.5 95.33
0.25 200 2.8 124.6
0.25 220 3.1 162.6
0.25 260 3.6 220.8
0.25 290 4.0 290.6
0.25 330 4.6 373.3
0.25 360 5.0 479.3
0.25 430 6.0 572.2
0.25 500 6.9 795.6
0.25 580 8.1 1018
0.25 700 9.7 1224
0.35 84 1.5 33.43
0.35 110 2.0 74.23
0.35 140 2.5 106.8
0.35 170 3.0 152.5
0.35 200 3.6 222.6
0.35 230 4.1 293.0
0.35 280 5.0 434.6
0.35 340 6.1 672.0
0.35 400 7.1 863.7
0.35 450 8.0 1053
Table3.3:ResultsofEkin,Emag,andfdipforMHDdynamosatri/ro=0.15undertheFTboundarycondition. ri/roRa[×103 ]Ra/RacritEkinEmagΛRmfdipfmagfitM[ZAm2 ] 0.154404.0297.13.2440.0324121.9−−− 0.157607.0823.75.1930.0519202.9−−− 0.158708.0801.4501.65.02200.20.4941.43570.06 0.159809.0846.4579.05.79205.70.5161.86675.20 0.15110010.1748.4444.74.45193.40.3490.86027.69 0.15130011.91493300.53.01273.20.1170.32213.04 0.15150013.81867135.61.36305.50.1550.3919.155 0.15170015.62141234.32.34327.20.1720.42014.67
Table3.4:ResultsofEkin,Emag,andfdipforMHDdynamosatri/ro=0.25undertheFTboundarycondition. ri/roRa[×103 ]Ra/RacritEkinEmagΛRmfdipfmagfitM[ZAm2 ] 0.251401.976.691.491×10−4 1.491×10−6 61.68−−− 0.251602.259.05958.69.5954.330.8602.116270.4 0.251802.562.89109710.9756.080.8672.397286.4 0.252002.885.49844.48.4465.380.7841.828221.9 0.252203.1103.9769.27.6972.080.7572.441200.0 0.252603.6236.448.40.484108.70.6442.47734.87 0.252904.0289.383.00.830120.30.6202.61041.30 0.253304.6248.6753.87.54111.50.6022.287103.3 0.253605.0297.7769.07.69122.00.5621.935146.6 0.254306.0455.6491.04.91150.90.5221.88783.72 0.255006.9642.2305.23.05179.20.4561.55162.10 0.255808.1906.8126.01.26212.90.4121.55650.84 0.257009.7122125.40.254247.1−−−
Table3.5:ResultsofEkin,Emag,andfdipforMHDdynamosatri/ro=0.35undertheFTboundarycondition. ri/roRa[×103 ]Ra/RacritEkinEmagΛRmfdipfmagfitM[ZAm2 ] 0.35841.535.092.376×10−3 2.376×10−5 41.89−−− 0.351102.043.61819.68.2046.700.8164.713297.6 0.351402.589.35140814.166.830.7243.174356.0 0.351703.0106.8950.29.5073.080.7394.239236.8 0.352003.6136.3890.28.9082.550.6923.900200.3 0.352304.1193.1938.49.3898.260.6322.946203.4 0.352805.0311.9895.68.96124.90.5562.848162.3 0.353406.1508.7713.67.14159.50.4802.006132.7 0.354007.1837.673.460.735204.70.3762.07124.92 0.354508.0106010.610.106230.2−−−
Fig. 3.1:The kinetic energy density, as a function of the Rayleigh number, calculated in spherical shells with different geometries. Red, green, and blue lines indicate the linear fitting forri/ro = 0.15,0.25,and0.35, respectively.
Fig. 3.2:Time evolution of energy density in various amplitudes of the initial magnetic field atRa/Racrit= 2.5atri/ro = 0.25.
Fig. 3.3:Time evolution of energy density in the case of sustained dynamo at Ra/Racrit = 2.8 atri/ro = 0.25. The red and blue lines denote the kinetic and magnetic energy density, respectively.
Fig. 3.4:The kinetic and magnetic energy density as a function of the Rayleigh number in spherical shells with different geometries. The black, red, and blue points are theEkin values in the non-MHD cases,Ekin values in the MHD cases, and Emagvalues in the MHD cases, respectively. The ”F” denotes the failed dynamo cases.
Fig. 3.5:The Elsasser numberΛas a function of the Rayleigh number for different spher-ical shell geometries. The red, blue, and green circles indicateLambdavalues for the cases ofri/ro= 0.15,0.25,and0.35, respectively.
Fig. 3.6:The magnetic Reynolds numberRmas a function of the Rayleigh number for different spherical shell geometries. The red, blue, and green circles indicate Rmvalues for the cases ofri/ro = 0.15,0.25,and0.35, respectively.
Fig. 3.7:Spectrum of the kinetic and magnetic energy density for the case with Ra/Racrit = 2.0 at ri/ro = 0.35 and lmax = 47. The red and blue lines denote the kinetic and magnetic energy density, respectively.
Fig. 3.8:Spectrum of the kinetic and magnetic energy density for the case with Ra/Racrit = 2.0 at ri/ro = 0.35 and lmax = 95. The red and blue lines denote the kinetic and magnetic energy density, respectively.
Fig. 3.9:Spatial pattern of the flow and magnetic fields for the case withRa/Racrit = 11.9atri/ro = 0.15. The radial magnetic field,Br, at the CMB, viewed from ϕ = π/2and3π/2, are plotted in the upper left and right, respectively. Thez component of the vorticity,ωZ, and magnetic field,BZ, at the equatorial plane are plotted in the lower left and right, respectively.
Fig. 3.10:Spatial pattern of the flow and magnetic fields for the case withRa/Racrit = 3.1atri/ro = 0.25. The radial magnetic field,Br, at the CMB, viewed from ϕ=π/2and3π/2, are plotted in the upper left and right, respectively. Thez component of the vorticity,ωZ, and magnetic field,BZ, at the equatorial plane are plotted in the lower left and right, respectively.
Fig. 3.11:Spatial pattern of the flow and magnetic fields for the case withRa/Racrit = 3.0atri/ro = 0.35. The radial magnetic field,Br, at the CMB, viewed from ϕ=π/2and3π/2, are plotted in the upper left and right, respectively. Thez component of the vorticity,ωZ, and magnetic field,BZ, at the equatorial plane are plotted in the lower left and right, respectively.
Fig. 3.12:Spectrum of the kinetic and magnetic energy density att/τη = 2.0for the case with Ra/Racrit = 11.9 atri/ro = 0.15. The red and blue lines denote the kinetic and magnetic energy density, respectively.
Fig. 3.13:Spectrum of the kinetic and magnetic energy density att/τη = 2.0for the case with Ra/Racrit = 3.1 at ri/ro = 0.25. The red and blue lines denote the kinetic and magnetic energy density, respectively.
Fig. 3.14:Spectrum of the kinetic and magnetic energy density att/τη = 2.0for the case with Ra/Racrit = 3.0 at ri/ro = 0.35. The red and blue lines denote the kinetic and magnetic energy density, respectively.
Fig. 3.15:Dipolarity fdip as a function of the Rayleigh number for different spherical shell geometries. The red, blue, and green points indicate fmag fit values in the cases ofri/ro = 0.15,0.25,and0.35, respectively. The bar represents the standard deviation for calculatingfmag fitevery time step.
Fig. 3.16:Magnetic energy density at the CMB of simulation data in l = 1 and fitting value ofl= 1as a function of spherical harmonic degree forri/ro= 0.25and Ra/Racrit= 2.8.
Fig. 3.17:The ratio of the dipolar magnetic energy density at the CMB to the extrap-olation of l = 1, fmag fit, as a function of the Rayleigh number for dif-ferent geometries. The red, blue, and green points indicate the cases of ri/ro = 0.15,0.25,and 0.35, respectively. The bar represents the standard deviation.
Fig. 3.18:Magnetic energy density at the CMB of simulation data at l = 1 and fitting value ofl= 1as a function of spherical harmonic odd degree. (a) Dipolar case atRa/Racrit= 7.1atri/ro = 0.35, and (b) Non-dipolar case atRa/Racrit= 10.1inri/ro = 0.15.
Fig. 3.19:Dynamo regime inri/ro = 0.15,0.25,and0.35. Red circles, blue triangles, green squares, and black crosses represent strong dipolar, weak dipolar, non-dipolar, and failed dynamo cases, respectively.
Chapter 4
Fixed heat flux
4.1 Result of thermal convection
Spherically symmetric heat equation in a non-convective static state is given by
∇2T = 1 r2
d dr
( r2dT
dr )
= 0. (4.1)
Eq. (4.1) indicates that the temperature gradient is proportional to the inverse square of r as dT /dr = A/r2, where A is a constant, and therefore the heat flow Q = 4πr2 ×dT /dr is constant. When the temperature is fixed at the ICB (T = Ti at r = ri) and the CMB (T = To at r = ro), the radial profile of the temperature is obtained by
T(r) = rori(Ti−To) ro −ri
1
r −riTi−roTo
ro−ri . (4.2)
Using the normalization of the temperature by the temperature difference between the ICB and the CMB and the radial distance by the thickness of a spherical shell, we rewrite the equation of the temperature gradient as
dT
dr =−rori
r2 . (4.3)
Then, the temperature gradients at the ICB and the CMB are indicated by βi =
−ri/ro andβ0 =−ri/r1o, respectively.
Here we consider the situation that the heat flux is constant at the boundaries of the spherical shell and is fixed in time. In order to estimate the critical Rayleigh number (Racrit) in fixed heat flux (FF) simulations, we carried out simulations of non-magnetic thermal convection using the governing equations of (2.6) without the Lorentz force term, (2.7), and (2.8). We assumed that the temperature gradient at the CMB was fixed asβ0 and the temperature gradients at the ICB were set as βi. In the FF simulations, the Rayleigh number based on the temperature gradient (the flux Rayleigh number RaF) at the CMB is used as an input parameter. The flux Rayleigh numberRaF is defined by
RaF = αTg0⟨−dTdr|r=ro⟩L4
(ri/ro)νκT . (4.4)
perature difference between the ICB and the CMB,∆T, obtained in the simulation results. Hereafter, we use the conventional Rayleigh numberRaunless otherwise noted.
Table 4.1 summarizes the simulation results of the averaged Ekin in the time in-terval fromt/τν =4.5 to 6, corresponding to a quasi-steady state. The simulation results indicate thatEkin is proportional toRa, as we found in the fixed tempera-ture simulations like Chap. 3.1. The critical Rayleigh numbers are estimated to be Racrit = 7.04×104,4.26×104,and3.01×104 inri/ro = 0.15,0.25,and0.35, respectively, by the same method as Al-Shamali et al. (2004). Racrit of FF cases are smaller than that of FT cases in all radius ratios. These simulation results show that the thermal convection occurs easily in FF cases than in FT cases.
We calculated the averaged kinetic energy density Ekin in the time interval from t/τν = 4.5 to 6. Table 4.2 shows Ekin of larger Ra cases than those in Table 4.1. To compare FT and FF cases, Ekin of FT cases are shown again in Table 4.3. Fig.4.1 shows the averaged kinetic energy density as a function of the flux Rayleigh number in FT and FF cases at ri/ro = 0.15,0.25, and 0.35. Red and blue symbols correspond to FT and FF cases, respectively. First, we found that Ekin of FF cases is larger than that of FT cases. Because the critical Rayleigh number of FF cases is smaller than that of FT cases, FF cases cause more intense convection than FT cases at the same flux Rayleigh number. Second, we found that the difference ofEkin between FT and FF cases becomes small at large flux
Rayleigh number atri/ro = 0.25and0.35. To understand the obtained difference in detail, we calculated the radial profiles of the temperature gradient in FT and FF cases atri/ro = 0.25and0.35. Fig.4.2 shows the radial profiles of the absolute value of the temperature gradient at small (RaF = 2.0 × 105) and large flux Rayleigh numbers (RaF = 1.1×106) in FT and FF cases at ri/ro = 0.25,and 0.35, respectively. Blue and red symbols respectively correspond to the results of small and large RaF cases, and the dotted and solid lines indicate FT and FF cases, respectively. In the larger range ofRaF, the temperature gradient becomes similar in FT and FF cases at the both ri/ro, indicating that the difference of the boundary condition is not significant in the largeRaF range.