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!
Juliane

Hi, Juliane
Thank you for your reply.
I have changed the ./source/heating_model/radioactive_decay.cc 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
  end
end

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

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


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
  end
end

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
  end
end

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

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




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
end
# 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)
  end
end

subsection Formulation
  set Formulation = Boussinesq approximation
end


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
  end
  subsection Latent heat melt
    set Melting entropy change                =  -300
  end
  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
  end
  subsection Compositional heating
	set Compositional heating values          = 0, 1e-5, 0, 0, 0
  end 
end

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
  end

end


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
  end
end

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
  end
end
```![1635228866(1)|690x260](upload://AaorMiAUmzPmuR4lS8iNEj1IsFh.png)

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,
Juliane