The Non‐Negligible Effect of Viscosity Diffusion on the Geodynamo Process

Exploring geodynamo's scaling laws is of great importance when a tremendous gap still lies between realistic physics and estimated model parameters. Existing scaling laws need to be tested by numerical simulations. To boost these, taking the outer core viscosity ν , we studied its impacts on weak and strong field geodynamo outputs by varying ν in two orders of magnitudes. In the weak field mode, the fluid velocity u varies with ν by a scaling law of u∼ν0.52 . While in the strong field mode, u varies very slowly by a scaling law of u∼ν0.10 . The magnetic field B does not change much with ν when the driving force is not too strong ( Ra=12,000 ) but decreases with ν by a scaling law of B∼ν−0.20 when geodynamo operates in a very vigorous mode ( Ra=50,000 ). The reason that u increases with ν is essential that increasing ν breaks the Taylor‐Proudman constraint and drops the critical Rayleigh number, but this kind of increase of u does not give a stronger B . Furthermore, we conducted a local force balance analysis and demonstrated that the balance shifts under different ν . In a quasi‐MAC regime, the Lorentz force instead of the inertial force enters into the first‐order agostrophic force balance, though the viscous force still plays a role. Comparing with other scaling laws previous wisdom holds, we propose that the effect of viscosity diffusion is non‐negligible for both existing weak and strong dynamo regimes. Still, appropriate assumptions can be made to take this variation into account.


Introduction
Paleomagnetic records indicate that the Earth's geomagnetic field has existed for at least 3 billion years, with a roughly reversal rate of hundreds of thousands of years. Earth's magnetic field is generated and maintained by a strongly coupled mechanism known as geodynamo (Roberts, 2007). As emerging from numerical simulations, details of these mechanisms achieve a self-sustaining magnetic field (Aubert, 2020;Buffett, 2000;Yadav et al., 2015). From 1950 to 1990, the dynamo theory has shown that the principle of a homogeneous dynamo process is possible and has elucidated the parameter prerequisites and model constraints. In 1995, there was a breakthrough in the numerical simulation of magnetohydrodynamic (MHD) flow by Glatzmaiers and Roberts, called GR95 (Glatzmaiers & Roberts, 1995b). The GR95 model showed that such models could describe and explain many basic properties of the geomagnetic field, including dipole-dominated morphology, reversals, and westward drift (Glatzmaiers & Roberts, 1995a). Furthermore, they studied viscous and electromagnetic coupling between the inner core and the outer core, which brought differential rotation of the inner core (Glatzmaiers & Roberts, 1995a). A couple of years later, numerous numerical geodynamo models emerged, for a comprehensive review see Christensen and Wicht (2007). For those models, the magnetic field properties closely match the geomagnetic field in terms of spatial spectra and magnetic field morphology, secular variation, and occasionally the characteristics of dipole reversals (Aubert, 2020;Calkins et al., 2017;Wicht & Sanchez, 2019;Yadav et al., 2016).
Nevertheless, the huge gap between the simulations' parameter and their expected values in the Earth's core still exists (Aubert, 2018;Schaeffer et al., 2017;Schwaiger et al., 2019). Most models have four critical dimensionless parameters: the Rayleigh number a R , measures thermal driven force; the Ekman number E, measures the relative importance of viscous forces to Coriolis forces; the Prandtl number r P , measures the ratio of viscosity to thermal diffusivity; and the magnetic Prandtl number m P , measures the ratio of viscosity to magnetic diffusivity . Among the unsolved issues, the value space of controlling Abstract Exploring geodynamo's scaling laws is of great importance when a tremendous gap still lies between realistic physics and estimated model parameters. Existing scaling laws need to be tested by numerical simulations. To boost these, taking the outer core viscosity  , we studied its impacts on weak and strong field geodynamo outputs by varying  in two orders of magnitudes. In the weak field mode, the fluid velocity u varies with  by a scaling law of   0.52 u . While in the strong field mode, u varies very slowly by a scaling law of   0.10 u . The magnetic field B does not change much with  when the driving force is not too strong (  12,000 a R ) but decreases with  by a scaling law of    0.20 B when geodynamo operates in a very vigorous mode (  50,000 a R ). The reason that u increases with  is essential that increasing  breaks the Taylor-Proudman constraint and drops the critical Rayleigh number, but this kind of increase of u does not give a stronger B. Furthermore, we conducted a local force balance analysis and demonstrated that the balance shifts under different  . In a quasi-MAC regime, the Lorentz force instead of the inertial force enters into the first-order agostrophic force balance, though the viscous force still plays a role. Comparing with other scaling laws previous wisdom holds, we propose that the effect of viscosity diffusion is non-negligible for both existing weak and strong dynamo regimes. Still, appropriate assumptions can be made to take this variation into account.

DONG ET AL.
So, the fundamental purpose of second-generation geodynamo models is to explore the parameter space and derived scaling laws (Aubert, 2020;Christensen & Aubert, 2006;Schaeffer et al., 2017). Previous studies have shown promising results with regards to this problem, for example, by analyzing the effects of one single parameter: a R (Aubert et al., 2017;Davidson, 2013;Kuang et al., 2008;Olson et al., 1999;Sreenivasan et al., 2014;Wang et al., 2013), E (Kuang et al., 2017;Sarson et al., 1998), and m P (Calkins et al., 2017;Christensen & Aubert, 2006;Schaeffer et al., 2017;Simitev & Busse, 2005), or by performing comprehensive study in parameter space (Christensen & Aubert, 2006;Christensen et al., 1999;Davidson, 2013;Gillet & Jones, 2006;Jones, 2011). Many scaling laws and quantitative conclusions have thus been derived. Jones (2011) investigated in which parameter space does the geodynamo behaves as dipole-dominated or multipole-dominated. Starchenko and Jones (2002) discovered how to estimate the strength of the flow and the magnetic field from a given parameter. Christensen and Aubert (2006) studied whether the magnetic field's power was related to various diffusion coefficients. However, the four major dimensionless parameters mentioned above usually depend on all of the physical coefficients for different geodynamo models. Thus, their variations are not self-dependable, which leads to the doubt on the applicability of existing conclusions. Besides, some scaling laws related to additional dimensionless parameters which make the scaling rules more physically complex. It is significantly difficult to understand the physical mechanisms behind them.
A majority of works about Earth's core's viscosity focus on the effects on the geodynamo process. Sarson et al. (1998) studied the impact of the variation of E on the dynamo and found that both the magnetic and the flow field are concentrated near the ICB. Soderlund et al. (2012) highlighted the significant role of viscosity in their numerical simulation. Magnetic forces did not seem to play a major role, contrary to what is expected for the real Earth. King and Buffett (2013) confirmed the non-negligible of viscosity diffusion in their simulation. Cheng and Aurnou (2016) showed that the diffusionless scaling laws were hiding an actual dependency upon viscosity. Oruba and Dormy (2014) derived alternative scaling statutes in which the magnetic field intensity depends upon viscosity and rotation rate. In conclusion, viscosity plays a significant role in the geodynamo process. However, previous studies do not directly determine the quantitative impact of viscosity, nor does it provide an intuitive physical mechanism.
To investigate the effects of viscosity on the geodynamo mechanism and deduce the output when the estimated numerical model parameters are close to the physical parameters, we adopt the Modular, Scalable, Self-consistent, and Three-dimensional (MoSST) model and conduct numerical experiments to address these issues. MoSST model has been developed from the KB97 model (Kuang & Bloxham, 1997). In addition to some improvements in the numerical schemes, it also includes some additional functional modules/features. In most of the numerical geodynamo models in practice, the viscous diffusion time is introduced as a typical timescale. While in the MoSST model, the magnetic diffusion time is used. Of all the dimensionless parameters, only the Ekman number contains the viscous term. Therefore, the advantage of using MoSST model to study the influence of viscosity is that changing E can be directly equivalent to change  , without DONG ET AL. This study is structured as follows: The MoSST model is introduced in the Section 2. In the Section 3, we study the effects of viscosity on fields and the length scale in the outer core and analyze the system's force balance. In the Section 4, the physical mechanisms behind the geodynamo model are investigated. Finally, the summary and discussions are delivered in Section 5.

The MoSST Model
The geodynamo model's governing equations in the spherical coordinates are composed of the momentum equation (Navier-Stokes equation), the magnetic induction equation, the energy equation (conduction/ diffusion equation of thermal field), and the nondivergent constraints of both the velocity and magnetic fields. Under the Boussinesq approximation, dimensionless forms of the above equations are (Kuang & Bloxham, 1997) where o R is the magnetic Rossby number, u is the flow velocity in the outer core, P is the pressure after the correction, B is the magnetic induction intensity, E is the Ekman number, a R is the Rayleigh number, Θ is the temperature perturbation, r is the radius,  q is the modified Prandtl number, and 0 T is the conducting background state (   0 u B ).
The dimensionless parameters in Equation 1 are defined as: which  is the magnetic diffusion, Ω is the angular velocity, o r is the radius of the outer core,  is the kinematic viscosity,  is the thermal expansion coefficient in the outer core, T h is the temperature gradient at ICB, g is gravitational acceleration, and  is the thermal diffusivity.
It is worthwhile to mention that, of the four parameters in Equation 2,  appears only in the definition of E, which is significantly different from the dimensionless description of other models (Kono & Roberts, 2002). Therefore, this dynamo parameter definition allows easy isolation of the viscous effect by changing only the Ekman number. Meanwhile, taking magnetic diffusion time     2 / o r as the typical time scale is direct and appropriate. Other nondimensional scales are: (Kuang & Bloxham, 1999 where [] denotes the difference across the boundaries, E is the dimensionless electrical field. A D″ layer with 1 / 10 conductivity of the CMB is included outside the CMB. Fixed heat flux boundary condition     / 0, r at both ICB and CMB, is adopted for the temperature field. The initial state of our model is derived from previous simulation results.
The MoSST model is discretized by the fourth-order compact finite difference scheme in the radial direction, using zero points of Chebyshev polynomial expansion as the radial mesh configuration nodes. The number of nodes is 35, 39 and 19 for the inner core, outer core, and D″ layer, respectively. On the spherical DONG ET AL. . According to the truncation orders, the number of spherical direction nodes is 50 at  direction and 64 at  direction. The total number of meshes in the simulation is 277,830.
To avoid high truncation for the initial transient period and to reduce the CPU time for the process, hyperdiffusion for viscosity  is introduced by (Kuang & Bloxham, 1999): where  0 is the original viscosity,  is an artificially selected small number, l is the spherical harmonic order, and 0 l is an artificially selected spherical harmonic order. Similarly, thermal hyperdiffusion and magnetic hyperdiffusion are also introduced. This approach was first introduced by Glatzmaiers and Roberts (1995a), who obtained the first three-dimensional self-consistent dynamo model of GR95. The hyperdiffusion scheme of Equation 3 is a relatively weaker form than GR95. We choose   0 4, 0.05 l  in our simulation. The impacts of hyperdiffusion on our scaling laws will be discussed in Section 4.4.

Numerical Simulation Results
Dimensionless parameters control numerical simulation results. The MoSST model has four dimensionless parameters:  Table 1. The outer core values estimated in Table 1 are calculated from the physical parameters estimated by Olson (2007) and Christensen and Wicht (2007). The four parameters indicate a large gap between the numerical model and the realistic Earth. Unfortunately, a high-performance parallel machine's current computing power is far below the capacity required to reach Earth parameter values. So far, significant efforts have to be made in the future.

Diagnostic Outputs
Since the governing equations are strongly nonlinear systems with highly coupled physical fields, huge iteration steps are needed to reach a stable electromagnetism output state after the initial state is input. To determine whether the system has been self-evolved to this state, we rely on the diagnostic outputs. An example of the diagnostic output plot is shown in Figure 1, with the horizontal coordinate as the magnetic free decay time, measured by the free diffusion time as    The ratio of thermal diffusivity to magnetic diffusivity 6 10 1 *Leaving a R blank here because the super adiabatic heat flow (or, more generally, the buoyancy flux) is not well constrained for the Earth's core yet.

Table 1
Definitions of Dimensionless Parameters in MoSST perturbation (Θ), respectively. The typical spatial scales are corresponding to the three physical fields, velocity field length scale ( u l ), magnetic field length scale ( B l ), and temperature perturbation length scale ( Θ l ).
With these definitions, we estimate the stable output basing on the fact that all physical quantities and their typical scales have been stabilized for at least one  f . The final result is an average over the stable period, as indicated by the red line in Figure 1. The transient time depends on the parameter selection, boundary conditions, initial state, and hyper diffusion scheme. As limited by the computing capacity, in general, it takes 4-5 days to obtain a stable output. For some extreme cases, it may need 2 weeks or more to receive one sound output.

Variations of Physical Fields With Ekman Number
By changing the dimensionless parameter E when  12,000 a R , a series of stable output results are obtained, as shown in Table 2. According to the former study of Kuang and Bloxham (1999)  and more close to the real value in Earth's core ( 9 10 ), m a R , S a R and the following force balance may be different with those of Kuang and Bloxham (1999). Since the geodynamo system operates at a very nonlinear state, and the initial condition may affect the final results, we have selected the suitable outputs and get rid of the singular data set when searching for and analyzing the regularity. As shown in Table 2, the flow velocity u increases significantly with E, and its amplitude varies by more than one order. Similarly, the temperature perturbation Θ increases monotonically with E, but the change is very small ( 6%), suggesting that even a tiny temperature perturbation can lead to a significant increase in flow velocity. By contrast, the magnetic field B changes nonmonotonically with E, and it appears a trend of up first (from The results also show that a substantial increase in flow velocity does not necessarily result in considerable magnetic field growth.
In Table 2, we also list the typical spatial scales of the physical fields. As can be seen, u l increases monotonically with E by a variation of 53%. For the typical length scale of B, B l does not change monotonically with E. B l does not change much, and its variation is as small as  8%. In contrast to the Θ variation, the temperature perturbation length scale Θ l decreases monotonically with E, and the amplitude is less than 16%.
In our model, we take  3,500 km o r and  1,200 km i r , so the flow field's typical scale is equivalent to 2-3 wavelengths in the outer core. From our numerical results, there are ∼10 wavelengths of the magnetic field in the outer core, and the scale is ∼1/3-1/5 of the flow field. The temperature perturbation is about half the wavelengths in the outer core, larger than the convection's typical scale. The specific scales of these three physical fields are all shown to be typically in large scales. The spatial scale of the temperature perturbation is larger than the scale of convection, and the scale of convection is larger than the scale of the magnetic field. These may reflect the fact that the temperature perturbation motivates convection, and convection motivates the magnetic field. Details of these dynamics at different scales, and the stacking of energy at various scales, can be found in the work of Huguet and Amit (2012) and Calkins et al. (2015).  Table 2. It is clear that the influence of viscosity on the geodynamo occurs mainly in the flow field: the higher the viscosity, the more vigorous the convection. For quantitative studies, former scaling relations fitted the numerical results with power laws, such as   f E (Christensen & Aubert, 2006;King & Buffett, 2013). In this form, both fields and their typical length scales would vanish in the inviscid limit (  0 E ). Thus, here we assume a different scaling law of the form: where 0 f is the independent of E, and represents its inviscid limit. Then the least square fit is used to estimate the three parameters 0 f , , and .
By applying Equation 4 to the velocity field u and its typical length scale u l , we have

Force Balance
The fluid motion in the outer core is affected by various forces. And the force equilibrium mechanisms determine the dynamics in the outer core.
In this part, we analyze details of the force balance under different viscosities. The momentum equation in Equation 1 gives Equation 8 includes five forces: inertial force I F , Coriolis force C F , pressures gradient P F , Lorentz force L F , viscous force V F , and thermal buoyancy A F (Archimedes force).
Among them, the leading-order force balance is between C F and P F , with the remaining puny part compensated by others, which is called the quasi-geostrophic (QG) state (Schwaiger et al., 2019). If L F is close to C F and far beyond other forces ( , ,and A I V F F F ), this will be called the magneto-quasi-geostrophic state, in that the uncompensated part of  C P F F (also called the ageostrophic Coriolis) is mainly balanced by L F at zeroth order. Traditionally, geodynamo is assumed to operate on a QG-MAC state (shorted for MAC), in which the ageostrophic C F is balanced by  A L F F at first order. A MAC balance is usually treated as equal to the socalled strong field dynamo, while a discrepancy exists between the two definitions. In strong field dynamos, the ratio between L F and C F defined as the Elsasser number (in dimensional form) is is of   , B is averaged as 0.36 for all Ekmans, then we have  Λ 0.13, which supports that they are weak field dynamos, but the possibility that a MAC balance exists cannot be ruled out according to Λ. The force balance according to Λ could be misleading (Soderlund et al., 2012) since that is only a rough order estimation. The real force balance should be based on each force's exact calculation in Equation 8 with the numerical results. By following Kuang (1999), we use the local force defined by (taking L F as example)  always about three orders smaller than the Coriolis force and at least one order smaller than other forces for all E, which proves that the systems are all far from MAC state. Furthermore, we divide the force balance into four regions: Ⅰ, Ⅱ, Ⅲ, and Ⅳ. In region Ⅰ, the ageostrophic C F is mainly compensated by I F (cyan) and A F , which is called the QG-CIA balance (shorted for CIA). In region Ⅱ, with the fast increase of I F , the ageostrophic C F is mainly compensated by I F , which is called the QG-CI balance (shorted for CI). Similarly, we call region Ⅲ the CIV and region Ⅳ the CV balance. Though fluctuating, V F (blue) is generally increasing and becomes more and more significant in force balance with increasing of E, which is consistent with the increase of the fluid core's viscosity. When E decreases toward the real value for physical Earth, the system is expected to enter into a CIA balance in this weak field regime.
To check if there is a local distribution of the MAC state, we calculated the local effective Elsasser number  Λ and effective Rayleigh number  R defined as (Kuang, 1999): The curls in the above two equations are used to get the ageostrophic force. When   Λ 0.5 means that the ageostrophic Lorentz force is comparative to the ageostrophic Coriolis force, and it is similar for  R. Snapshot of  Λ and  R distribution for  12,000  . Figure 4b shows the distribution of  Λ with grid points at different r, , and Figure 4c shows the histogram distribution of  Λ. whose axis coincides with the rotation axis and tangents with the inner core at ICB). Both of  Λ and  R show symmetry about the equator to a large extent. Figure 4b shows the distribution of  Λ with grid points of different r (from io r to o r ),  (from 0 to  / 2). There are some relatively obvious  Λ near the ICB. The histogram distribution of  Λ is plotted in Figure 4c, which shows that most areas have a typical  Λ of about   3 2 10 . Figure 4 indicates that the ageostrophic Coriolis force is mainly compensated by A F outside TC but by I F inside TC. L F only plays a negligible role in the force balance in this weak field regime.
The ratio between magnetic energy density ( m E ) and kinematic energy density ( k E ) is another criterion to evaluate dynamo states. In a strong field dynamo, usually one has  m k E E (at least two orders). The energy ratio is calculated by We use the numerical results of Table 2 to calculate it and use Equation 4 to fit its variation, then get Figure 5. As can be seen, the energy ratio almost decreases monotonically with E by a fitted exponent of 0.94, which results from increasing of u but not decreasing of B (see Figure 2). Another noticed character is that the magnetic energy is always smaller than the kinematic energy, proving that the system is in a weak field regime and far from a MAC balance.

Field Configuration
To figure out the effect of viscosity variation on the dynamo system, we mapped the flow fields with respect to different E, as shown in Figure 6. Each figure contains two halves. The left-colored half represents the axisymmetric (  0 m ) toroidal flow, and the right half represents the axisymmetric poloidal flow. The results are all time-averaged outputs over one  f The toroidal flow ( T u ) and poloidal flow ( P u ) are defined as The left part of Figure 6 is the axisymmetric zonal component of velocity  u , the right part is the axisymmetric streamline (Ψ u ), defined by

a) c) b)
where S u is the axisymmetric velocity,    sin s r is the distance from the field point to Earth's rotation axis.
As can be seen from Figure 6, the mean flow velocity increases with E. When    6 2.5 10 E , the liquid iron flow mainly concentrates on the TC's profile. As E increases gradually, the flow concentrates on the TC and extends toward CMB on the equator plane. These two features are both reflected in the toroidal flow (left half). For the poloidal flow (right half), a significant increase of flow scales can be found both inside and outside the TC.
The axisymmetric magnetic field is shown in Figure 7. As can be seen, for the toroidal field, it could be locally strong, especially inside the TC, though the system is in a weak field regime as a whole. The toroidal field strength decreases and tends to concentrate in the TC with increasing of E. There appears no significant difference for the poloidal magnetic field, and the dipolar field outside CMB (where only a poloidal field exists) all dominate.

Our
 12,000 a R results are all proven to be in the weak field regime. To study the strong field mode, we gradually increase a R and two typical values,  30,000 a R and 50,000 are selected and listed in The exponential factor is about 0.10 for both equations, which is much smaller than that for the weak field regime    6 to 4), which indicates that the system may be in a strong field regime. Second, we calculated the RMS values of the local forces defined by Equation 10 and normalized all of them by C F . The results are shown in Figure 9. As can be seen, rvc increases monotonically with E and the other three force ratio fluctuate. With increasing of E, rvc varies from ∼15% to 50%. rmc and rvc are always comparable to each other, and they are about 20%. ric is as small as 5%, which indicates that inertial plays a tiny role in the force balance. According to the above variances of force ratios, Figure 9 where the ageostrophic C F is mainly compensated by   A L V F F F and Ⅲ-VC, where the ageostrophic C F is mainly compensated by V F . As E decreases toward the real value, the system is expected to be more close to a MAC state, which is consistent with geodynamo theoretical consensus.
We call the two solutions for  50,000 a R (    5 2.5 10 E and   5 5 10 ) in a quasi-MAC state, since V F still plays a role but will turn into more and more insignificant with decreasing of E.
Third, the local quasi-MAC balance is checked by the effective Elsasser number  Λ and effective Rayleigh number  R defined in Equations 11 and 12. Results are shown in Figure 10. As can be seen, the ageostrophic C F is mainly compensated by L F near ICB, by A F near CMB, and by V F near the rotational axis (refer to Figure 9). Figure 10b shows that large  Λ distributes at ICB and some part of CMB, while Figure 10c shows that  Λ has a typical value of ∼0.2.
Lastly, we check the energy density ration / m k E E , as defined by Equation 13, using the RMS u and B listed in Table 3. The result is shown in Figure 11. As can be seen, when E increases from   5 2.5 10 to   4 1.25 10 , the energy ratio decreases from 23.9 to 7.4. As E varies, the magnetic energy is larger than the kinematic energy by about an order. This again proves that the system is in a quasi-MAC balance. If E decreases to the real value, this ratio is expected to become even larger, consistent with MAC's theoretical predictions or strong field dynamo regime (Schaeffer et al., 2017). DONG ET AL. . Figure 10b shows the distribution of  Λ with grid points at different r, , and Figure 10c shows the histogram distribution of  Λ.

a) b)
c) Figure 11. The ratio between magnetic energy density and kinematic energy density for  50, 000 a R . Red and black lines are for numerical and fitted results, respectively.

Discussions
The outer core flow increases with viscosity, especially under a weak field, significantly different from the hydrodynamics system's general concept. Previous studies on geodynamo proposed the same conclusion. Here, we discuss the similarities and differences with previous studies and analyze the variation mechanism.

Influence of Viscosity on the Critical Driving Force and Flow Scales
The fluid motion in the outer core is strongly affected by Earth's rotation, and the rotation, therefore, induces a fundamental change in the dynamo pattern of the outer core. The leading-order momentum equation is the geostrophic equilibrium (which can also be seen in Figure 3, where Coriolis forces dominate throughout), where only Coriolis forces and the pressure gradient are retained for Equation 8 Taking the curl of both sides of Equation 20, we have This shows that dynamo flow is essentially two-dimensional and is constant along the rotation axis, namely the Taylor-Proudman (T-P) theorem. Taylor (1963) first proposed that the convection pattern is a series of parallel cylindrical convection rolls parallel to the rotation axis and rotating along the rotation axis. However, in a closed spherical system, this form of flow can neither transfer heat to the outside nor maintain the dynamo process. Therefore, the T-P constraint must be broken so that the governing equation's first-order term is ageostrophic. To study the influence of viscosity alone, by neglecting A F , L F and I F (e.g., Region Ⅳ in Figure 3 and Region Ⅲ in Figure 9), taking the curl of Equation 8 and obtaining the vorticity equation, and taking the radial component of the vorticity equation, then we obtain where    ω u is vorticity. From Equation 22, it is apparent that, under the  0 r u boundary conditions, larger E will result in larger r u . This may partly and qualitatively explain why u increase with E or  . u is also affected by other velocity components, together with the vorticity in the right-hand side of the equation and other forces. Therefore, u does not exhibit a simple proportional relationship with E. From another point of view, increasing the viscosity help break the T-P constraint by lowering the critical Rayleigh number a c R , Which marks the critical driving force for thermal convection in a fast-rotating system. Chandrasekhar (1981) introduced the steady state rotational thermal convection theory and presented that Some studies have also analyzed the convective roll diameter conv l , and found that r r is the spherical shell thickness Chandrasekhar, 1981;Johnson & Constable, 1998;Roberts & King, 2013). Figure  Although researches on E through a c R and conv l are many, only King and Buffett (2013) provided a quantitative result on the effect of E on u by  1/2 1/3 u C E . However, this result cannot be simply interpreted as the direct relationship between velocity and viscosity in the outer core because their work is based on is the Reynolds number and    3 3 0 C PH / A is the convection energy. If we assume that P and A are independent of  , we can derive:    1/6 u . This inference contradicts to the dynamo system. Therefore, the most significant contribution of this study is that it provides the only direct and reliable quantitative relationship between u and  to date by fitting numerical dynamo data. The variation of with  is weaker than other studies, especially under a weak field (exponential factor  0.12). The convective column results from the ideal approximation, and realistic dynamo convection is complex and affected by various forces. It can also be seen from the distribution of r u on the equatorial plane, as shown in Figure 12, where convection rolls are elongated and twisted in the radial direction. They are no longer an ideal cylindrical shape, and their spatial scale could not be simply represented by the diameter of convection rolls anymore.
When the T-P constraint is broken, the convection becomes quasi-three-dimensional. There is a radial component, whose amplitude is only 2% of  u (referred to Figure 6). Stress-free boundaries lead to larger zonal flows, so it might also explain why our results differ from others researches in that most studies use the no-slip boundaries conditions, which will be affected by the Ekman layer (Wicht & Sanchez, 2019;Yadav et al., 2016). However, this 2% radial portion is critical to the generation and maintenance of the geodynamo process. Were it not for this small r u , the heat cannot be transported outside, and the flow will degenerate DONG ET AL. into two-dimensional. According to Cowling's (1933) antidynamo principle, two-dimensional flow cannot maintain the dynamo process.

The Role of Viscosity on the Magnetic Field and Flow Field
Next, we focus on the determining factors for the magnitude of the magnetic field. As is mentioned, geodynamo is speculated to operate at a strong field regime, such that B satisfies    Λ 1 O . In this situation, the force balance determines the magnitude of B. Starchenko and Jones (2002) proposed that it is rather the driving force and force balance than diffusion coefficients that determine the flow velocity and the magnetic field. Christensen and Aubert (2006) assumed that the diffusion effects do not play a major role in the governing equations, and they found that the available energy (represented by the flow Rayleigh number a Q R ), rather than the force balance, determined the magnitude of the magnetic field. Christensen (2010) summarized scaling laws for the magnetic field and proposed that the magnitude of B depended on the available energy flux, which could overcome the ohmic loss. They argued that rotation only determines whether a dipole field is generated. Jones (2011)  , which partly supports the conclusion that the magnitude of B is not strongly associated with viscosity. While when, the magnetic field is annihilated, indicating that  has a significant impact on magnetic field generation upon a critical point (via c a R to m a R ).    7 6 10 E For quasi-MAC regime, the magnetic field decreases significantly with (decreases by 35.4% as E is 5 times larger), and the increase of flow is also noticeable (16.5%), indicating that the influence of viscosity on the flow and magnetic field cannot be ignored.

The Role of Viscosity on Geodynamo Force Balance
Geodynamo is expected to run in the MAC state. The existing geodynamo simulations can explain a couple of observations in principle. Those model parameters adopted in simulations are still far from physical. It leads to active arguments on whether existing geodynamo models can run in a good regime representing the physical Earth. After checking the exact force balance or fitting the numerical data, studies showed that viscosity, instead of the magnetic field, plays a major role in the force balance for most of the existing geodynamo simulation results (King & Buffett, 2013;Oruba & Dormy, 2014;Soderlund et al., 2012). Most existing dynamo models are, in fact running in VAC, not MAC balance (Yadav et al., 2016), and a real MAC state is hard to reach except using a more physical parameter (e.g.,  5,000 c a a R R,   7 10 E ) and simulating at very high resolution (Aubert, 2019;Schaeffer et al., 2017).
Variations of  can bring substantial changes to the force balance. In the weak field regime, as is shown in Figure 3, as E increases, the system gradually transfers from CIA mode to VC mode. CIA mode may result in a strongly driven but nonmagnetic system or multiple dynamo regimes (Soderlund et al., 2012). When    6 3 10 E , V F grows more and more obvious in the force balance. Though fluctuates with E, L F is always fairly weak (Figure 3), rendering it only plays a minor role in the force balance, and the dynamo state is far from a MAC balance (Figure 4). The energy ratio between m E and k E increases as E decreases, while we always have  m k E E in our parameter space ( Figure 5). In the quasi-MAC regime, as E increases, the system transfers from a quasi-MAC mode to VC mode (Figure 9). When    5 5 10 E , C F grows more and more important in the force balance. When    5 5 10 E , the system enters into a quasi-MAC regime, verified by the force ratio distribution (Figure 9), the local effective  R and  Λ distribution (Figure 10), and the energy ratio ( Figure 11). As E decreases toward the real value, L F plays a more and more important role in core dynamics, and the system is expected to evolve into a real MAC regime, as is estimated theoretically.
The exponent of scaling relation of u and  decreases from 0.52 for  12,000 a R (Figure 2) to 0.10 for  50,000 a R (Figure 8), which indicates that the effects of viscosity on fluid velocity is much smaller for strong-driven and in the quasi-MAC regime than the weakly driven and in the weak field regime. The DONG ET AL.

10.1029/2020JB021281
15 of 19 reason is not that the relative important of V F is smaller in a quasi-MAC balance (compare Figure 3 with Figure 9), but that L F enters into the first-order ageostophic balance ( Figure 10) and replaces I F , and much of the kinematic energy is transferred into the magnetic energy (compare Figure 5 with Figure 11).

The Effect of the Hyperdiffusion Assumption
Hyperdiffusion could be avoided, while the temporal and spatial resolution must be improved. To improve the computational efficiency, especially simulate at high a R and low E, hyperdiffusion is retained in this study. Hyperdiffusion mainly suppresses the development of small-scale convection and has no significant effect on large-scale flow, which is one of the major reasons that the typical spatial scale of various physical fields is all large in our simulations. Researches on the influences of hyperdiffusion are many. For example, Zhang and Jones (1997) found that hyperviscosity significantly increased the equivalent Ekman number and changed the convection kinetics. Grote et al. (2000) concluded that hyperdiffusion could affect the magnetic field generation mechanism and even reverse polarity.
In our simulations, we choose   To obtain quantitative properties and provide quantified "large-scale patterns" which will be used to derive the scaling laws, we should find the highest spherical harmonic degree H L , above which the hyperdiffusive effect will be dominant. Only use the part of the solutions with spherical harmonic coefficients of degrees  H l L in our analysis. With this approach, we would hope to find asymptotic behaviors of the "large-scale" dynamo solutions. Then we can see the effective Ekman number H E to account for the hyperdiffusivity used in the simulation, and derive the scaling laws with respect to H E . The relative importance of the hyperdiffusivity for large-scale solutions (  H l L ) can be described by the ratios of the following viscous dissipation terms, The H L and H E are calculated and listed in Table 4. As can be seen, in the weak field regime, the effective Ekman numbers are overall  10 30 times the Ekman numbers used in the analyses as mentioned earlier. Though this will affect the quantitative results for each E, it would not change the whole regularity (e.g., the exponential in the scaling law). The effects of hyperdiffusion are equivalent to multiply each E in our scaling laws by a factor of ∼20. The maximum degree for hyperdiffusion H L is about 5-7, and scales smaller than these (  H l L ) will be affected substantially. These effects should be kept in mind for references to this study.

Conclusions
To systematically study the effects of viscosity on the geomagnetic field dynamo both in a weak and strong field, we adopt the newly developed MoSST model to perform numerical experiments on the dynamo process. The advantage of the MoSST model is that, since the magnetic diffusion time is used as the typical time scale,  is only included in E, making the variation of E equivalent to that of  . Major conclusions reached through our numerical investigations are as follows 1. The role of core viscosity on dynamo is non-negligible and different for weak field and quasi-MAC regime. In a weak field regime, variations of  may affect the fluid velocity and its typical length scale  (Figure 7). The magnitude of B is not so sensitive to  . While in the quasi-MAC regime, the effect of  on u and u l is much smaller, but B decreases fast with  ( Figure 8). The force balance varies solely with the relative importance of V F (Figure 9). 2. As  decreases toward the real value, the dynamo system involves into CIA balance for the weak field regime and MAC balance for the quasi-MAC regime. Our scaling relations are different from existing scaling laws not only in the scaling form but also in the exponent. For example, in a CIA state, former studies (King & Buffett, 2013)  3. The drop of fitted exponent from 0.52 for the weak field to 0.10 for the quasi-MAC state is not attributed to the relative importance of V F becomes weaker (comparing Figures 3 and 9), but because that L F instead of I F plays a more important role in the latter regime.
There are limitations in this study. When changing the value of E and . a R The other two-dimensionless parameters ( o R and  q ) are kept constant. Changes in them may quantitatively affect the conclusions drawn here. Solutions of the quasi-MAC state maybe not sufficient, and the resolutions need to be improved. Hyperdiffusion still has effects on core dynamics especially for smaller length scales (  7 l ). We will aim to solve these in future research.