Questions about subsection " Heating model ": redioactive decay and latent heat melt

Hello, as the title,I have two questions about "Heating model "
(1) How to set three or more different radioactive decay in three or more compositions. For example, I want to set compositon 1 、2、3 with different radioactive decay concentrations, but I found “heating model/redioactive decay” options only have crust and mantle radioactive concentrations to set.
(2)I want to set a model : the melt does not move relative to the solid, but the melt have the latent heat release and density change.This couldn’t be achieved by postprocessor. Should I select the “latent heat melt” both in subsections " Heating model " and " Material "? And the first composition should be set to porosity or pyroxenite? I don’t understand the description about this setting.

I’d be grateful to you for answering my question.

Hi Qyunian,

For your first question, I would recommend using the “compositional heating” plugin. It implements a model in which the magnitude of internal heat production is determined from fixed values assigned to each compositional field.

For your second question, yes, you can use the “latent heat melt” material model, together with the “latent heat” heating model (the “latent heat melt” heating model assumes that you have a separate compositional field for the melt that can move relative to the solid, which it seems is not what you want in your model). The “latent heat melt” material model assumes that the material is peridotite if you do not have any compositional fields in the model. If compositional fields are specified, the model assumes that the first compositional field is the fraction of pyroxenite and the rest of the material is peridotite. So it does not matter how you call this compositional field, you just need to be aware of this convention.

I hope that helps!

Hi, Juliane
Thank you for your reply.
I have changed the ./source/heating_model/ and ./include/heating_model/radioactive_decay.h files to achieve my porpuse. But I found another question is about the output value of “radioactive_decay”. I thought the unit of this value is W/m3 (i.e. density*concentrations[element]*heating_rates[element]), but when I set multi-composition fields with different densitis, the output value won’t change for different density fields. I found this case will happen even I use the original radioactive_decay files(.cc&.h) , but when I change the background density , the output values will change but won’t show different values for different density fields.
Is this a bug in ASPECT?

# test average (density and radioactive decay value)

set Dimension                              = 2
set Maximum time step                      = 1e6
set Output directory                       = output-test-new
set Use years in output instead of seconds = true
set End time                               = 1.4e7

subsection Solver parameters
  subsection Stokes solver parameters
    set Linear solver tolerance = 1e-8
    set Number of cheap Stokes solver steps = 100

subsection Boundary velocity model
  set Tangential velocity boundary indicators = left, right, top, bottom

subsection Initial temperature model
  set List of model names = function
  subsection Function
    set Function expression     = 1000

subsection Boundary temperature model
  set Fixed temperature boundary indicators = left, right, top, bottom
  set List of model names = constant
  subsection Constant
    set Boundary indicator to temperature mappings  = left:1000, right:1000, bottom:1000, top:1000

subsection Geometry model
  set Model name = box
  subsection Box
    set X repetitions = 200
    set Y repetitions = 50
    set X extent      = 400
    set Y extent      = 100

subsection Gravity model
  set Model name = vertical
  subsection Vertical
    set Magnitude = 9.81

subsection Mesh refinement
  set Initial adaptive refinement        = 0
  set Initial global refinement          = 0
  set Time steps between mesh refinement = 0

subsection Compositional fields
  set Compositional field methods = field
  set Number of fields = 4
  set Names of fields = Density_1 , Density_2 , Radio_1, Radio_2
# porosity 

subsection Initial composition model
  set Model name = function
  subsection Function
    set Coordinate system   = cartesian
    set Variable names      = x , y 
    set Function expression = if(x>0 && x<200  && y<40, 1, 0) ; if(x>200 && x<400 && y<40 , 1, 0) ; if(x>0 && x<200  && y>60, 1 ,0) ; if(x>200 && x<400 && y>60, 1 ,0)

subsection Formulation
  set Formulation = Boussinesq approximation

subsection Heating model
  set List of model names =  radioactive decay
  subsection Function
    set Function expression                   = 4e-9 ; # 0 ; 4e-9          # without mantle, PKT, IBC_remain, MIC ; W/kg
  subsection Latent heat melt
    set Melting entropy change                =  -300
  subsection Radioactive decay
    set Crust composition number      		  = 2                          # radioactive decay for MIC
    set Crust defined by composition          = true
    set Heating rates                         = 0.1, 0.2		   # U, Th
    set Half decay times                      = 4.47e9, 14.00e9
    set Initial concentrations crust          = 100,  500          #
    set Initial concentrations mantle         = 10, 50
    set Number of elements                    = 2
  subsection Compositional heating
	set Compositional heating values          = 0, 1e-5, 0, 0, 0

subsection Material model

  set Model name =  diffusion dislocation

  subsection Diffusion dislocation
    set Reference temperature  = 1000
    set Reference viscosity    = 5e20 
    set Minimum viscosity      = 5e16
    set Maximum viscosity      = 5e26
    set Viscosity averaging scheme =  harmonic
    set Effective viscosity coefficient           = 1.0

    set Densities              = 10000, 20000, 30000, 40000, 40000                                 # mantle, PKT, IBC_remain, MIC_rich, MIC
    set Thermal diffusivity    = 1e-6                                                    
    set Thermal expansivities  = 3e-5
    set Heat capacity          = 1250

    set Prefactors for diffusion creep            = 1e-21          #e=0
    #set Prefactors for diffusion creep            = 2.09325e-18     #e=100
    #set Prefactors for diffusion creep            = 4.381709e-15
    #set Prefactors for diffusion creep            = 9.172027e-12
    set Prefactors for dislocation creep          = 1e-100
    set Activation energies for diffusion creep   = 0    #100kj
    set Activation energies for dislocation creep = 0
    set Activation volumes for diffusion creep    = 0
    set Activation volumes for dislocation creep  = 0
    set Stress exponents for diffusion creep      = 1
    set Stress exponents for dislocation creep    = 1

    set Grain size                                = 1      #defult:1e-3
    set Grain size exponents for diffusion creep  = 0      #defult:3.
    set Maximum strain rate ratio iterations      = 40
    set Minimum strain rate                       = 1.4e-20
    set Strain rate residual tolerance            = 1e-22


subsection Postprocess
  set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization, composition statistics
  subsection Visualization
    set List of output variables = density, viscosity, melt fraction, temperature anomaly, heating
    set Output format = vtu
    set Time between graphical output = 3e7

subsection Solver parameters
  set Temperature solver tolerance = 1e-12
  subsection Stokes solver parameters
    set Linear solver tolerance = 1e-8
    set Number of cheap Stokes solver steps = 0
    #set Stokes solver type                  = block GMG
    # A higher restart length makes the solver more robust for large viscosity contrasts.
    set Maximum number of expensive Stokes solver steps = 10000
    set GMRES solver restart length         = 1000

Dear Qyunian,

You are using the Boussinesq Approximation in your model. The Boussinesq Approximation assumes that the density is constant (and equal to the reference density) everywhere except in the buoyancy term. That means in the heat equation, the reference density is used (rather than the real density). This is probably why you do not see a change in the heat production rate when you change densities of individual compositions, but the heat production changes when you change the reference density.

Best regards,