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
orDissNb
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 variablesstrat
,polind
,g0
,g1
andg2
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
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 ofpenaltyfac
.
Transport properties¶
difExp (default
difExp=-0.5
) is a real. This is the exponent that is used whennVarVisc=2
,nVarDiff=2
ornVarCond=4
.
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.
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 whennVarCond=1
.con_FuncWidth (default
con_FuncWidth=0.25
) is a real. This parameter is only used whennVarCond=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:Spherical harmonic degree \(\ell\)
Spherical harmonic order \(m\)
Real amplitude (\(\cos\) contribution)
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 ass_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 inn_impS
(the maximum allowed is 20), and it’s given by the number ofsCMB
defined in the same namelist. The default value ofimpS
is zero (no entropy disturbance). If it is set in the namelist for an integer greater than zero, thensCMB
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 ofs_bot
ands_top
). Each four consecutive numbers are interpreted as follows:Highest amplitude value of the entropy boundary condition, stored in array
peakS(20)
. WhenimpS<0
,peakS
is a relative amplitude in comparison to the \((\ell=0,m=0)\) contribution (for example, the cases_top= 0 0 -1 0
).\(\theta\) coordinate (input has to be given in degrees), stored in array
thetaS(20)
.\(\phi\) coordinate (input has to be given in degrees), stored in array
phiS(20)
.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:Spherical harmonic degree \(\ell\)
Spherical harmonic order \(m\)
Real amplitude (\(\cos\) contribution)
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 asxi_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 inn_impXi
(the maximum allowed is 20), and it’s given by the number ofxiCMB
defined in the same namelist. The default value ofimpXi
is zero (no chemical composiiton disturbance). If it is set in the namelist for an integer greater than zero, thenxiCMB
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 ofxi_bot
andxi_top
). Each four consecutive numbers are interpreted as follows:Highest amplitude value of the chemical composition boundary condition, stored in the array
peakXi(20)
. WhenimpXi<0
,peakXi
is a relative amplitude in comparison to the \((\ell=0,m=0)\) contribution (for example, the casexi_top= 0 0 -1 0
).\(\theta\) coordinate (input has to be given in degrees), stored in array
thetaXi(20)
.\(\phi\) coordinate (input has to be given in degrees), stored in array
phiXi(20)
.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.