Physical parameters namelist

This namelist contains all the appropriate relevant control physical parameters.

Dimensionless control parameters

  • ra (default ra=0.0) is a real. This the thermal Rayleigh number expressed by

    \[Ra = \frac{\alpha g_o \Delta T d^3}{\kappa\nu}\]
  • raxi (default raxi=0.0) is a real. This the compositional Rayleigh number expressed by

    \[Ra_\xi = \frac{\alpha g_o \Delta \xi d^3}{\kappa_\xi\nu}\]
  • ek (default ek=1e-3) is a real. This is the Ekman number expressed by

    \[E = \frac{\nu}{\Omega d^2}\]
  • pr (default pr=1.0) is a real. This is the Prandtl number expressed by

    \[Pr = \frac{\nu}{\kappa}\]
  • sc (default sc=10.0) is a real. This is the Schmidt number expressed by

    \[Sc = \frac{\nu}{\kappa_\xi}\]
  • prmag (default prmag=5.0) is a real. This is the magnetic Prandtl number expressed by

    \[Pm = \frac{\nu}{\lambda}\]
  • po (default po=0.0) is a real. This is the Poincaré number expressed by

    \[Po = \frac{\Omega_p}{\Omega}\]
  • prec_angle (default prec_angle=23.5) is a real. This is the angle between the precession and the rotation axes expressed in degrees.

  • radratio (default radratio=0.35) is a real. This is the ratio of the inner core radius \(r_i\) to the outer core radius \(r_o\):

    \[\eta = \frac{r_i}{r_o}\]
  • strat (default strat=0.0) is a real. This is the number of density scale heights of the reference state:

    \[N_\rho = \ln \frac{\tilde{\rho}(r_i)}{\tilde{\rho}(r_o)}\]
  • DissNb (default DissNb=0.0) is a real. This is the dissipation number:

    \[Di = \frac{\alpha_o g_o d}{c_p}\]

    Warning

    This can only be provided as a replacement input of strat. I.E., when one wants to define a reference state, one has to specify either strat or DissNb in the input namelist.

  • polind (default polind=1.5) is a real. This is the polytropic index, which relates the background temperature to the background density:

    \[\tilde{\rho} = \tilde{T}^m\]

    Warning

    Be careful: in its current version the code only handles adiabatic backgrounds, therefore changing polind physically means that the nature of the fluid (in particular its Grüneisen parameter) will change. For an ideal gas, it actually always follows \(m+1=\frac{\gamma -1}{\gamma}\)

  • l_isothermal (default l_isothermal=.false.) is a logical. When set to .true., makes the temperature background isothermal (i.e. \(\tilde{T}=cst.\)). In that case, the dissipation number \(Di\) vanishes and there is no viscous and Ohmic heating left. The only difference with the Boussinesq set of equations are thus restricted to the density background \(\tilde{\rho}\) and its radial derivatives that enters the viscous stress. This approximation is also called the zero Grüneisen parameter and was extensively explored by Denise Tortorella during her PhD.

Heat sources and sinks

  • epsc0 (default epsc0=0.0) is a real. This is the volumetric heat source \(\epsilon_0\) that enters the thermal equilibrium relation:

    (1)\[-\nabla\cdot\left(\tilde{\rho}\tilde{T}\nabla s\right) + \epsilon_0\,f(r)=0\]

    The radial function \(f(r)\) can be modified with the variable nVarEps that enters the same input namelist.

  • epscxi0 (default epscxi0=0.0) is a real. This is the volumetric source \(\epsilon_\xi\) that enters the compositional equilibrium relation:

    (2)\[-\nabla\cdot\left(\tilde{\rho}\nabla \xi\right) + \epsilon_\xi=0\]
  • nVarEps (default nVarEps=0) is an integer. This is used to modify the radial-dependence ofthe volumetric heat source, i.e. \(f(r)\) that enters equation (1).

    nVarEps=0

    Constant, i.e. \(f(r)=\hbox{cst.}\).

    nVarEps=1

    Proportional to density, i.e. \(f(r)=\tilde{\rho}(r)\).

    nVarEps=2

    Proportional to density times temperature, i.e. \(f(r)=\tilde{\rho}(r)\tilde{T}\).

Realistic interior models

  • interior_model (default interior_model="None") is a character string. This defines a polynomial fit of the density profile of the interior structure of several astrophysical objects. Possible options are "earth", "jupiter", "saturn" and "sun" (the naming is not case sensitive).

    Warning

    When interior_model is defined the variables strat, polind, g0, g1 and g2 are not interpreted.

    The subroutine radial gives the exact details of the implementation.

  • r_cut_model (default r_cut_model=0.98) is a real. This defines the cut-off radius of the reference model, i.e. the fluid domain is restricted to radii with \(r\leq r_{cut}\).

The following input parameters will thus define a polynomial fit to the expected interior structure of Jupiter until 99% of Jupiter’s radius (assumed here at the 1 bar level)

interior_model="JUP",
r_cut_model   =0.99e0,

Gravity

The radial dependence of the gravity profile can be adjusted following

(3)\[g(r)=g_0+g_1\,\frac{r}{r_o}+g_2\,\left(\frac{r_o}{r}\right)^2\]

The three following parameters are used to set this profile

  • g0 (default g0=0) is the pre-factor of the constant part of the gravity profile, i.e. \(g_0\) in equation (3).

  • g1 (default g1=1) is the pre-factor of the linear part of the gravity profile, i.e. \(g_1\) in equation (3).

  • g2 (default g2=0) is the pre-factor of the \(1/r^2\) part of the gravity profile, i.e. \(g_2\) in equation (3).

Centrifugal acceleration

The centrifugal acceleration can be computed for a polytropic background

  • dilution_fac (default dilution_fac=0.0) is the ratio of the centrifugal acceleration at the equator to the surface gravitational acceleration.

    (4)\[m=\frac{\Omega^2 d}{g_o}\]

Phase field

  • stef (default stef=1.0 is a real. This is the Stefan number (used when the phase field is plugged in). It is expressed by the ratio of the latent heat per unit mass associated with the solid-liquid transition and the specific heat:

    \[St = \frac{\mathcal{L}}{c_p\Delta T}\]
  • tmelt (default tmelt=0.0) is a real. This is the dimensionless melting temperature.

  • epsPhase (default epsPhase=0.01) is a real. This is the dimensionless interface thickness between the solid and the liquid phase (sometimes known as the Cahn number).

  • phaseDiffFac (default phaseDiffFac=1.0) is a real. This is a coefficient that goes in front of the diffusion term in the phase field equation.

  • penaltyFac (default penaltyFac=1.0) is a real. This is coefficient used for the penalisation of the velocity field in the solid phase. The smaller the coefficient, the stronger the penalisation. Since this is a nonlinear term, it is handled explicitly and the time step size should be decreased with the square of penaltyFac.

Transport properties

Electrical conductivity

There are several electrical conductivity profiles implemented in the code that can be chosen with the nVarCond input variable. The following one corresponds to a constant electrical conductivity in the deep interior (\(r<r_m\)) and an exponential decay in the outer layer.

(5)\[\begin{split}\sigma(r)=1+ (\sigma_m-1)\left(\frac{r-r_i}{r_m-r_i}\right)^a \quad \hbox{for}\quad r<r_m, \\ \sigma(r)=\sigma_m \exp \left[a \left(\frac{r-r_m}{r_m-r_i}\right)\frac{\sigma_m-1}{\sigma_m}\right] \quad\hbox{for}\quad r\geq r_m.\end{split}\]
  • nVarCond (default nVarCond=0) is an integer. This is used to modify the radial-dependence of the electrical conductivity.

    nVarCond=0

    Constant electrical conductivity, i.e. \(\sigma=\hbox{cst.}\)

    nVarCond=1

    \(\sigma\propto\tanh[a(r-r_m)]\)

    nVarCond=2

    See equation (5).

    nVarCond=3

    Magnetic diffusivity proportional to \(1/\tilde{\rho}\), i.e.
    \[\lambda=\frac{\tilde{\rho}_i}{\tilde{\rho}}\]

    nVarCond=4

    Radial profile of the form:
    \[\lambda=\left(\frac{\tilde{\rho}(r)} {\tilde{\rho}_i}\right)^{\alpha}\]
  • con_RadRatio (default con_RadRatio=0.75) is a real. This defines the transition radius \(r_m\) that enters equation (5).

  • con_DecRate (default con_DecRate=9) is an integer. This defines the decay rate \(a\) that enters equation (5).

  • con_LambdaMatch (default con_LambdaMatch=0.6) is a real. This is the value of the conductivity at the transition point \(\sigma_m\) that enters equation (5).

  • con_LambdaOut (default con_LambdaOut=0.1) is a real. This is the value of the conduvity at the outer boundary. This parameter is only used when nVarCond=1.

  • con_FuncWidth (default con_FuncWidth=0.25) is a real. This parameter is only used when nVarCond=1.

  • r_LCR (default r_LCR=2.0) is a real. r_LCR possibly defines a low-conductivity region for \(r\geq r_{LCR}\), in which the electrical conductivity vanishes, i.e. \(\lambda=0\).

Thermal diffusivity

  • nVarDiff (default nVarDiff=0) is an integer. This is used to change the radial-dependence of the thermal diffusivity:

    nVarDiff=0

    Constant thermal diffusivity \(\kappa\)

    nVarDiff=1

    Constant thermal conductivity, i.e.
    \[\kappa =\frac{\tilde{\rho}_i}{\tilde{\rho}(r)}\]

    nVarDiff=2

    Radial profile of the form:
    \[\kappa=\left(\frac{\tilde{\rho}(r)} {\tilde{\rho}_i}\right)^{\alpha}\]

    nVarDiff=3

    polynomial-fit to an interior model of Jupiter

    nVarDiff=4

    polynomial-fit to an interior model of the Earth liquid core

Viscosity

  • nVarVisc (default nVarVisc=0) is an integer. This is used to change the radial-dependence of the viscosity:

    nVarVisc=0

    Constant kinematic viscosity \(\nu\)

    nVarVisc=1

    Constant dynamic viscosity, i.e.
    \[\nu =\frac{\tilde{\rho}_o}{\tilde{\rho}(r)}\]

    nVarVisc=2

    Radial profile of the form:
    \[\nu=\left(\frac{\tilde{\rho}(r)} {\tilde{\rho}_i}\right)^{\alpha}\]

    where \(\alpha\) is an exponent set by the namelist input variable difExp.

Anelastic liquid equations

Warning

This part is still work in progress. The input parameters here are likely to be changed in the future.

  • epsS (default epsS=0.0) is a real. It controls the deviation to the adiabat. It can be related to the small parameter \(\epsilon\):

    \[\epsilon \simeq \frac{\Delta T}{T} \simeq \frac{\Delta s}{c_p}\]
  • cmbHflux (default cmbHflux=0.0) is a real. This is the CMB heat flux that enters the calculation of the reference state of the liquid core of the Earth, when the anelastic liquid approximation is employed.

  • slopeStrat (default slopeStrat=20.0) is a real. This parameter controls the transition between the convective layer and the stably-stratified layer below the CMB.

Boundary conditions

Thermal boundary conditions

  • ktops (default ktops=1) is an integer to specify the outer boundary entropy (or temperature) boundary condition:

    ktops=1

    Fixed temperature (Boussinesq) or entropy (anelastic) at outer boundary: \(s(r_o)=s_{top}\)

    ktops=2

    Fixed temperature gradient (Boussinesq) or entropy gradient at outer boundary: \(\partial s(r_o)/\partial r = q_t\)

    ktops=3

    Only use it in anelastic models: fixed temperature at outer boundary: \(T(r_o)=T_{top}\)

    ktops=4

    Only use it in anelastic models: fixed temperature gradient at outer boundary: \(\partial T(r_o)/\partial r = q_t\)

  • kbots (default ktops=1) is an integer to specify the inner boundary entropy (or temperature) boundary condition.

  • s_top (default s_top= 0 0 0.0 0.0) is a real array of lateraly varying outer heat boundary conditions. Each four consecutive numbers are interpreted as follows:

    1. Spherical harmonic degree \(\ell\)

    2. Spherical harmonic order \(m\)

    3. Real amplitude (\(\cos\) contribution)

    4. Imaginary amplitude (\(\sin\) contribution)

    For example, if the boundary condition should be a combination of an \((\ell=1,m=0)\) sherical harmonic with the amplitude 1 and an \((\ell=2,m=1)\) spherical harmonic with the amplitude (0.5,0.5) the respective namelist entry could read:

    s_top = 1, 0, 1.0, 0.0, 2, 1, 0.5, 0.5, ! The comas could be left away.
    
  • s_bot (default s_bot=0 0 0.0 0.0) is a real array. This is the same as s_top but for the bottom boundary.

  • impS (default impS=0) is an integer. This is a flag to indicate if there is a localized entropy disturbance, imposed at the CMB. The number of these input boundary conditions is stored in n_impS (the maximum allowed is 20), and it’s given by the number of sCMB defined in the same namelist. The default value of impS is zero (no entropy disturbance). If it is set in the namelist for an integer greater than zero, then sCMB has to be also defined in the namelist, as shown below.

  • sCMB (default sCMB=0.0 0.0 0.0 0.0) is a real array of CMB heat boundary conditions (similar to the case of s_bot and s_top). Each four consecutive numbers are interpreted as follows:

    1. Highest amplitude value of the entropy boundary condition, stored in array peakS(20). When impS<0, peakS is a relative amplitude in comparison to the \((\ell=0,m=0)\) contribution (for example, the case s_top= 0 0 -1 0).

    2. \(\theta\) coordinate (input has to be given in degrees), stored in array thetaS(20).

    3. \(\phi\) coordinate (input has to be given in degrees), stored in array phiS(20).

    4. Angular width (input has to be given in degrees), stored in array widthS(20).

Boundary conditions for chemical composition

  • ktopxi (default ktopxi=1) is an integer to specify the outer boundary chemical composition boundary condition:

    ktopxi=1

    Fixed composition at outer boundary: \(\xi(r_o)=\xi_{top}\)

    ktopxi=2

    Fixed composition gradient at outer boundary: \(\partial \xi(r_o)/\partial r = q_t\)

  • kbotxi (default ktopxi=1) is an integer to specify the inner boundary chemical composition boundary condition.

  • xi_top (default xi_top= 0 0 0.0 0.0) is a real array of lateraly varying outer chemical composition boundary conditions. Each four consecutive numbers are interpreted as follows:

    1. Spherical harmonic degree \(\ell\)

    2. Spherical harmonic order \(m\)

    3. Real amplitude (\(\cos\) contribution)

    4. Imaginary amplitude (\(\sin\) contribution)

    For example, if the boundary condition should be a combination of an \((\ell=1,m=0)\) sherical harmonic with the amplitude 1 and an \((\ell=2,m=1)\) spherical harmonic with the amplitude (0.5,0.5) the respective namelist entry could read:

    xi_top = 1, 0, 1.0, 0.0, 2, 1, 0.5, 0.5, ! The comas could be left away.
    
  • xi_bot (default xi_bot=0 0 0.0 0.0) is a real array. This is the same as xi_top but for the bottom boundary.

  • impXi (default impXi=0) is an integer. This is a flag to indicate if there is a localized chemical composition disturbance, imposed at the CMB. The number of these input boundary conditions is stored in n_impXi (the maximum allowed is 20), and it’s given by the number of xiCMB defined in the same namelist. The default value of impXi is zero (no chemical composiiton disturbance). If it is set in the namelist for an integer greater than zero, then xiCMB has to be also defined in the namelist, as shown below.

  • xiCMB (default xiCMB=0.0 0.0 0.0 0.0) is a real array of CMB chemical composition boundary conditions (similar to the case of xi_bot and xi_top). Each four consecutive numbers are interpreted as follows:

    1. Highest amplitude value of the chemical composition boundary condition, stored in the array peakXi(20). When impXi<0, peakXi is a relative amplitude in comparison to the \((\ell=0,m=0)\) contribution (for example, the case xi_top= 0 0 -1 0).

    2. \(\theta\) coordinate (input has to be given in degrees), stored in array thetaXi(20).

    3. \(\phi\) coordinate (input has to be given in degrees), stored in array phiXi(20).

    4. Angular width (input has to be given in degrees), stored in array widthXi(20).

Boundary conditions for phase field

  • ktopphi (default ktopphi=1) is an integer to specify the boundary condition of the phase field at the outer boundary.

    ktopphi=1

    Fixed phase field at outer boundary: \(\phi(r_o)=\phi_{top}\)

    ktopphi=2

    Fixed phase field gradient : \(\partial \phi(r_o)/\partial r = 0\)

  • kbotphi (default kbotphi=1) is an integer to specify the boundary condition of the phase field at the inner boundary.

Mechanical boundary conditions

  • ktopv (default ktopv=2) is an integer, which corresponds to the mechanical boundary condition for \(r=r_o\).

    ktopv=1

    Stress-free outer boundary for \(r=r_o\):
    \[\begin{split}W_{\ell m}(r=r_o)=0, \quad \frac{\partial}{\partial r}\left(\frac{1}{r^2\tilde{\rho}} \frac{\partial W_{\ell m}}{\partial r}\right)=0 \\ \frac{\partial}{\partial r}\left(\frac{1}{r^2\tilde{\rho}} Z_{\ell m}\right)=0\end{split}\]

    ktopv=2

    Rigid outer boundary for \(r=r_o\):
    \[\begin{split}W_{\ell m}=0,\quad \frac{\partial W_{\ell m}}{\partial r}=0, \\ Z_{\ell m}=0\end{split}\]
  • kbotv (default kbotv=2) is an integer, which corresponds to the mechanical boundary condition for \(r=r_i\).

Magnetic boundary conditions

  • ktopb (default ktopb=1) is an integer, which corresponds to the magnetic boundary condition for \(r=r_o\).

    ktopb=1

    Insulating outer boundary:
    \[\frac{\partial g_{\ell m}}{\partial r}+\frac{\ell}{r}\,g_{\ell m}=0,\quad \frac{\partial h_{\ell m}}{\partial r}=0\]

    ktopb=2

    Perfect condutor:
    \[g_{\ell m} = \frac{\partial^2 g_{\ell m}}{\partial r^2}=0,\quad \frac{\partial h_{\ell m}}{\partial r}=0\]

    ktopb=3

    Finitely conducting mantle

    ktopb=4

    Pseudo-vacuum outer boundary:
    \[\frac{\partial g_{\ell m}}{\partial r}=0,\quad h_{\ell m}=0\]
  • kbotb (default kbotb=1) is an integer, which corresponds to the magnetic boundary condition for \(r=r_i\).

    kbotb=1

    Insulating inner boundary:
    \[\frac{\partial g_{\ell m}}{\partial r}-\frac{\ell+1}{r}\,g_{\ell m}=0,\quad \frac{\partial h_{\ell m}}{\partial r}=0\]

    kbotb=2

    Perfectly-conducting inner core:
    \[g_{\ell m} = \frac{\partial^2 g_{\ell m}}{\partial r^2}=0,\quad \frac{\partial h_{\ell m}}{\partial r}=0\]

    kbotb=3

    Finitely conducting inner core

    kbotb=4

    Pseudo-vacuum outer boundary:
    \[\frac{\partial g_{\ell m}}{\partial r}=0,\quad h_{\ell m}=0\]

Boundary condition for spherically-symmetric pressure

  • ktopp (default ktopp=1) is an integer, which corresponds to the boundary condition for the spherically-symmetric pressure at \(r=r_o\).

    ktopp=1

    The integral of the spherically-symmetric density perturbation vanishes.

    ktopp=2

    The spherically-symmetric pressure fluctuation vanishes at the outer boundary.