I’m trying to run a model in Aspect that has radioactive decay, and also uses two compositional fields with different densities, but I am noticing some strange behavior.
I start with a two layer model, with a basal layer that has slightly higher density. I also use the radioactive decay heating model, where the heating elements are enriched (by 20x) in the basal layer.
However, it appears that the initial density distribution is used to compute the heating throughout the simulation, instead of computing the heating with the updated density distribution as the two layers mix.
Initially the simulation starts with two layers, and after some time it looks like this:
The heating elements (top plots) generally partition to follow the compositional fields (bottom plots). However, there is a strange jump in the heat production that corresponds to the initial location of the basal layer. The change in heating is not that large, I had to adjust the scaling to see it. It is the transition from light red to dark red across the bottom of the plume.
It is visible in the background mantle as well, when the scale bar is properly adjusted.
My suspicion is that while it is using the compositional fields to properly determine the W/kg heating rate for each compositional field, it is using the initial density distribution to convert W/kg to the W/m^3 heating rates shown in the plots.This would explain why the differences in heating rate are small – the density difference between the top/bottom layer is only 59 kg/m^3.
Instead, aspect should use the density distribution at each timestep to convert the heating rate.
Does this look like the proper diagnosis? Any ideas on how to fix this?
(More images to follow - I could only put one image per post as a new user)
And here is the simulation after some time, but with the scalebar adjusted to see the jump in heating elements in the background mantle (transition from dark blue to light blue), at the same location as the original density difference:
Hi Robert - Thank you for posting to the forum! When you have a chance, can you post the relevant sections of the parameter file (heating model, compositional fields) or the whole parameter file? It is helpful to have that as a reference when looking at the relevant source code.
The source code (radioactive_decay.cc) looks like it is using the correct density, but this heating model was setup to handle a specific number of compositional fields.
Thanks for the reply!. Here are the relevant sections of the parameter file:
subsection Heating model
set List of model names = radioactive decay
subsection Radioactive decay
set Crust composition number = 0
set Crust defined by composition = true
set Number of elements = 4
# U38, U35, Th32, K40
set Heating rates = 9.4946e-5, 5.6840e-4, 2.6368e-5, 2.8761e-5
set Half decay times = 4.468e9, 0.704e9, 14.00e9, 1.248e9
set Initial concentrations mantle = 17.0e-3, 5.0e-3, 30.0e-3, 160.0e-3
######################################################################
# CHANGE: PPM CONCENTRATIONS ON BASAL LAYER (CRUST)
set Initial concentrations crust = 240.0e-3, 72.0e-3, 850.0e-3, 3700.0e-3
######################################################################
end
# set List of model names = compositional heating
# subsection Compositional heating
# set Compositional heating values = 6.039e-9, 1.2078e-07
# set Use compositional field for heat production averaging = 1, 1
# end
end
subsection Compositional fields
set Number of fields = 1
set Names of fields = basal_layer
set Compositional field methods = particles
set Mapped particle properties = basal_layer:initial basal_layer
end
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
##########################################################
# CHANGE: THICKNESS OF SHEEP LAYER
set Function constants = Thickness_layer=300e3
##########################################################
set Function expression = if(y<Thickness_layer, 1, 0)
end
end
We also set the density difference in the material model:
subsection Material model
set Model name = depth dependent
subsection Depth dependent model
set Base model = simple
set Depth dependence method = List
set Depth list = 670e3,2890e3
set Viscosity list = 2.05e21,6.15e22
end
subsection Simple model
##########################################################
# CHANGE: VISCOSITY AMBIENT MANTLE
# CHANGE: CONTRAST BETWEEN MATERIALS
set Viscosity = 2.05e21
set Composition viscosity prefactor = 1.
##########################################################
set Thermal expansion coefficient = 1e-5
set Thermal viscosity exponent = 4.8 # same as Li et al. 2018
set Reference density = 3300
set Reference temperature = 1600
##########################################################
# CHANGE: REFERENCE SURFACE DENSITY FOR THE MANTLE
# CHANGE: DENSITY DIFFERENCE BETWEEN LAYERS
set Density differential for compositional field 1 = 59
##########################################################
set Thermal conductivity = 4.7
end
set Material averaging = harmonic average
end
You can also download the entire parameter file here.
When the constant heating model is used (commented out), everything works as expected.
Hi Robert,
this looks weird indeed, but I think it is related to your model setup. You use the boussinesq approximation, which uses a reference profile of density for the temperature equation. Usually this reference density would be constant, but you have a compositional field with a higher density in the bottom layer, so if you use the default option for the reference profile (compute profile) it will use the initial conditions to compute this reference profile. And since this initial condition includes the high density bottom layer you get a jump in the reference density for your temperature equation. You can fix this by prescribing instead of computing the reference profile. Add the following to your parameter file:
subsection Adiabatic conditions model
set Model name = function
subsection Function
set Function constants =
set Function expression = T;P;RHO
set Variable names = depth
end
end
Replace T,P, and RHO by the reference temperature (constant for BA), pressure (rho*g), and density (constant for BA). This should hopefully fix the issue.
Thanks Rene! It did work correctly using this, except I put rho*g*depth for the pressure.
I see now how we didn’t notice this when using the compositional heating model, because the heating rates are input as W/m^3 so it never references the density when computing the heating rate. But now I’m curious, is the other density dependent term in the temperature equation (rho*Cp*…) using the incorrect (compute profile) reference density as well? I.e., even when using the compositional heating model we should use the Adiabatic conditions model as you describe?
Thanks Rene! It did work correctly using this, except I put rho * g * depth for the pressure.
Ah, of course.
Uh, yes I think you are right, the temperature equation uses the reference density for all terms. I think when we made the decision to include the compositional fields in the computation of the reference profile we had cases in mind where the composition is averaged from two fields (e.g. a mixture of harzburgite and basalt), in which case ignoring compositions would create a wrong reference profile. Using the new setting should fix the problem, and I also saw that the compute profile plugin (the default) has a setting called set Composition reference profile = initial composition that you can chance to set Composition reference profile = function, which lets you define a function for the reference composition (in your case leaving this function as 0.0 - the default - would also fix the issue).
But in your case the default setting for the composition reference profile is clearly not the right option. Any idea how we make this more clear in the documentation?
Perhaps a note in the Formulation section (p300, A.40) describing how the Adiabatic conditions model needs to be changed when using the BA with varying compositional densities. Another note could be added to the entry on p29 (2.11.5). Although, it would be better if when the BA was used the code automatically switched the reference profile to use a single constant density, as defined in the material model. Or at least, if a warning is output if BA is being used with compositional densities with the wrong Adiabatic conditions model. But it is a specific case that might not be common enough to warrant that.
