Dear developers
In a future paper I absolutely need to cool my CMB using the dynamic core temperature model:
###################################################
# SUBSECTION: Boundary temperature model
###################################################
# The CMB temperature will evolve and cool
subsection Boundary temperature model
set Fixed temperature boundary indicators = inner, outer
set List of model names = dynamic core
###################################################
# SUBSECTION: Dynamic core
###################################################
subsection Dynamic core
set Inner temperature = 4412 # We explore 4 different initial temperatures: 5912, 5412, 4912 and 4412 (Nakagawa and Tackley; 2010)
set Outer temperature = 300 # We keep this temperature fixed at 300 K.
set Core density = 12.5e3
set Gravity acceleration = 10.682 # Gravity value at CMB according to the PREM model (ascii file used in ASPECT).
set CMB pressure = 0.14e12
set Initial light composition = 0.0526211 # wt for light elements (Zhang et al. 2016).
set Max iteration = 30000
set Core heat capacity = 840
set K0 = 4.111e11
set Alpha = 1.35e-5 # Thermal expansivity for the core
set Core conductivity = 64 # Experimental data from Pozzo et al. (2022).
subsection Geotherm parameters
set Tm0 = 1695
set Tm1 = 10.9
set Tm2 = -8.0
set Theta = 0.11
set Composition dependency = true
set Use BW11 = false
end
###################################################
# SUBSECTION: Radioactive heat source
###################################################
subsection Radioactive heat source
# K in core
set Number of radioactive heating elements = 1
set Heating rates = 2.8761e-5 # Same value used in Citron et al. 2020
set Half life times = 1.25 # Unit: Gyr
set Initial concentrations = 4.0e-4 # Unit: ppm - We try 0, 200, 400, 600, 800 ppm.
end
end
end
the aim should be that, while the CMB cool, the solidus condition for the inner core should be satisfied, but they are never…
In the test directory, values are:
subsection Dynamic core
set Inner temperature = 4100
set Outer temperature = 300
set Core density = 12.5e3
set Gravity acceleration = 9.81 # For the core
set CMB pressure = 0.14e12
set Initial light composition = 0.042 # wt% for light elements
set Max iteration = 30
set Core heat capacity = 840
set K0 = 4.111e11
set Alpha = 1.35e-5 # Thermal expansivity for the core
set Core conductivity = 50
subsection Geotherm parameters
set Tm0= 1695
set Tm1= 10.9e-12
set Tm2= -8.0e-24
set Theta= 0.11
set Composition dependency= false
end
They are in agreement with Nimmo et al. (2004), the paper used to do this .cc
https://academic.oup.com/gji/article/156/2/363/2046212
If I use these e-12 and e-24 values, i got an inner core, but in an anomalous way: it start to grow from the beginnig, which is wrong according to the literature, and the T_CMB drop very slowly.
If you check the manual there is written that:
Parameter name: Tm0
Default value: 1695.
Pattern: [Double 0…MAX_DOUBLE (inclusive)]
Documentation: Melting curve ([Nimmo et al., 2004] eq. (40)) parameter Tm0. Units: \si{\kelvin}.
Parameter name: Tm1
Default value: 10.9
Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]
Documentation: Melting curve ([Nimmo et al., 2004] eq. (40)) parameter Tm1. Units: \si{\per\tera\pascal}.
Parameter name: Tm2
Default value: -8.0
Pattern: [Double -MAX_DOUBLE…MAX_DOUBLE (inclusive)]
Documentation: Melting curve ([Nimmo et al., 2004] eq. (40)) parameter Tm2. Units: \si{\per\tera\pascal\squared}.
So values should be INTEGERS as input, and not like the test. dynamic_core.cc agrees (?)
prm.enter_subsection("Geotherm parameters");
{
prm.declare_entry ("Tm0","1695.",
Patterns::Double (0.),
"Melting curve (\\cite{NPB+04} eq. (40)) parameter Tm0. Units: \\si{\\kelvin}.");
prm.declare_entry ("Tm1","10.9",
Patterns::Double (),
"Melting curve (\\cite{NPB+04} eq. (40)) parameter Tm1. "
"Units: \\si{\\per\\tera\\pascal}.");
prm.declare_entry ("Tm2","-8.0",
Patterns::Double (),
"Melting curve (\\cite{NPB+04} eq. (40)) parameter Tm2. "
"Units: \\si{\\per\\tera\\pascal\\squared}.");
prm.declare_entry ("Theta","0.11",
Patterns::Double (),
"Melting curve (\\cite{NPB+04} eq. (40)) parameter Theta.");
prm.declare_entry ("Composition dependency","true",
Patterns::Bool (),
"If melting curve dependent on composition.");
prm.declare_entry ("Use BW11","false",
Patterns::Bool (),
"If using the Fe-FeS system solidus from Buono \\& Walker (2011) instead.");
}
but if I use integers, not only i CANNOT, NEVER get an inner core (even with a 300 K surface and T CMB temperature), the only thing that it’s been fixed is that the temperature drops fine (another simulation):
I can’t figure out where is my problem. The cooling and inner core growth should be similar to this:
https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2010GC003031
If i plot the equation of the paper, i’d get a full freezing of the core at around 3000K
I attach the updated dynamic_core.cc (which is stucked in pull request).
dynamic_core.zip (9.3 KB)
and here a minimum working example (just to check if the inner core instantaneously form or not):
# A simple setup for convection in a quarter of a 2d shell. See the
# manual for more information.
set Dimension = 2
set Use years in output instead of seconds = true
set End time = 1.5e9
set Output directory = output-shell_simple_2d
subsection Material model
set Model name = simple
set Material averaging = harmonic average only viscosity
subsection Simple model
set Thermal expansion coefficient = 4e-5
set Viscosity = 1e22
end
end
subsection Boundary velocity model
set Tangential velocity boundary indicators = bottom, top
end
###################################################
# SUBSECTION: Geometry model
###################################################
# We run the simulation in a quarter of a shell with periodic boundaries.
subsection Geometry model
set Model name = spherical shell
subsection Spherical shell
set Inner radius = 3481000
set Opening angle = 90
set Outer radius = 6371000
end
end
###################################################
# SUBSECTION: Boundary temperature model
###################################################
# The CMB temperature will evolve and cool
subsection Boundary temperature model
set Fixed temperature boundary indicators = inner, outer
set List of model names = dynamic core
###################################################
# SUBSECTION: Dynamic core
###################################################
subsection Dynamic core
set Inner temperature = 2000 # We explore 4 different initial temperatures: 5912, 5412, 4912 and 4412 (Nakagawa and Tackley; 2010)
set Outer temperature = 300 # We keep this temperature fixed at 300 K.
set Core density = 12.5e3
set Gravity acceleration = 10.682 # Gravity value at CMB according to the PREM model (ascii file used in ASPECT).
set CMB pressure = 0.14e12
set Initial light composition = 0.1 #0.01 #0.0526211 # wt for light elements (Zhang et al. 2016).
set Max iteration = 30000
set Core heat capacity = 840
set K0 = 4.111e11
set Alpha = 1.35e-5 # Thermal expansivity for the core
set Core conductivity = 64 # Experimental data from Pozzo et al. (2022).
subsection Geotherm parameters
set Tm0 = 1695
set Tm1 = 10.9
set Tm2 = -8.0
set Theta = 0.11
set Composition dependency = true
set Use BW11 = false
end
###################################################
# SUBSECTION: Radioactive heat source
###################################################
subsection Radioactive heat source
# K in core
set Number of radioactive heating elements = 1
set Heating rates = 2.8761e-5 # Same value used in Citron et al. 2020
set Half life times = 1.25 # Unit: Gyr
set Initial concentrations = 0 # Unit: ppm - We try 0, 200, 400, 600, 800 ppm.
end
end
end
subsection Heating model
set List of model names = shear heating
end
###################################################
# SUBSECTION: Initial temperature model
###################################################
# The temperature at the boundaries evolves according to the "half space cooling model" and we set an age of 5e8 years to give a temperature contrast to the simulation and initiate convection as soon as possible.
subsection Initial temperature model
set Model name = function
# It's given a 200 K harmonic perturbation to initiate convection.
subsection Function
set Coordinate system = spherical
set Variable names = r, theta
set Function constants = r_inner=3481000, r_outer=6371000, w=2*pi, p=200.0, frequency=25
set Function expression = p * sin(pi * (r - r_inner) / (r_outer - r_inner)) * sin(frequency * w * theta)
end
end
subsection Boundary temperature model
set Fixed temperature boundary indicators = 0,1
end
subsection Gravity model
set Model name = ascii data
end
subsection Mesh refinement
set Initial global refinement = 5
set Initial adaptive refinement = 4
set Strategy = temperature
set Time steps between mesh refinement = 15
end
subsection Postprocess
set List of postprocessors = visualization, velocity statistics, temperature statistics, heat flux statistics, depth average, core statistics
subsection Visualization
set Output format = vtu
set Time between graphical output = 1e6
set Number of grouped files = 0
end
subsection Depth average
set Time between graphical output = 1e6
end
end
subsection Solver parameters
subsection Stokes solver parameters
set Stokes solver type = block GMG
end
end
So… Where is the problem? In the code? I misunderstood something? Please let me know and thank you in advance.
