Mass conservation and latent heating during a phase transition

ASPECT
1

I’m running models of crustal eclogitization that include a dynamic phase transition. My model is based on the Visco plastic material model.

What is the best way to make sure my phase transition conserves mass due to the density increase? (e.g., https://academic.oup.com/gji/article/186/1/6/2104990?login=true)

To conserve mass, is it enough to have a Compositing material model where the compressibility is assigned to the “multicomponent compressible” model, and everything else to “visco plastic”?

Relatedly, how can I incorporate latent heating into this reaction? I think that ASPECT can calculate this based on Clapeyron slopes (aspect/latent_heat.cc at edcb0fa34bfdf4a5a49d58a90112837dfb95501e · geodynamics/aspect · GitHub)

On the user’s end, is it enough to include “latent heat” in my list of heating models (along with shear heating, etc.), or do I need to add it as a third material model under Compositing? It seems that all the necessary information is already contained in the Visco Plastic model.

Apologies if this is covered somewhere, but I couldn’t find it in the manual explicitly and was having difficulty tracing the code to see what does what.

2

Hi Mitchell,

As you point out, it’s not a simple problem, and no, it’s not enough to use the compressibility from the “multicomponent compressible” model. We discuss the whole problem in more detail here https://doi.org/10.1093/gji/ggaa078, specifically, we also describe an example concerning phase transitions. The parameter files we used can be found in this repo GitHub - gassmoeller/formulations-of-compressible-mantle-convection-data, and the one in the folder slab-dehydration-reaction is the one with the phase transition.

Specifically, the problem is that even if you use a compressible model, the mass conservation equation needs to compute the density gradient somehow. Most implemetations in compressible mantle convection codes assume that the density follows some sort of reference profile (and we also have this option: The Anelastic Liquid Approximation). However, that means that the mass only changes in the correct way if you are on this reference profile (in terms of temperature, composition etc.). So if you whole model is just basalt/eclogite, that might be okay, but if you have other compositions, this will not work. So what we do instead is that we have the option to interpolate the density onto a finite-element field and then compute the density gradient from that, which always takes into account the actual density in the model. To use that option, you need to pick

subsection Formulation
  set Mass conservation = projected density field
  set Temperature equation = real density
end

in your input file. You also need to have a compositional field with the name projected_density that is used for the purpose of interpolating the density, for example:

subsection Compositional fields
  set Number of fields = 2
  set Names of fields = crust, projected_density
  set Compositional field methods = field, prescribed field
end

But for this to work, the material model needs to actually provide the outputs that can be interpolated onto the finite element field. The visco-plastic model does not do that at the moment (since it was written to be incompressible). It would need to have a function like this:

template <int dim>
void
NameOfMaterialModel<dim>::create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
{
  if (out.template get_additional_output<PrescribedFieldOutputs<dim>>() == NULL)
    {
      const unsigned int n_points = out.n_evaluation_points();
      out.additional_outputs.push_back(
        std::unique_ptr<MaterialModel::AdditionalMaterialOutputs<dim>>
        (new MaterialModel::PrescribedFieldOutputs<dim> (n_points, this->n_compositional_fields())));
    }
}

that creates the outputs, and then fill them in the evaluate function:

  if (PrescribedFieldOutputs<dim> *prescribed_field_out = out.template get_additional_output<PrescribedFieldOutputs<dim>>())
    for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
      prescribed_field_out->prescribed_field_outputs[i][density_field_index] = out.densities[i];

For an exmaple, you can have a look at benchmarks/compressibility_formulations/plugins/compressibility_formumations.cc, and the steinberger material model has this as well. So you would have to add this to the material model you want to use.

For the latent heat, the material model needs to compute terms that describe how much the energy should change, and I do not think that the viscoplastic/multicomponent compressible model computes these (but I could be wrong, so maybe someone who actually uses these models can correct me). The latent heat material model is one that computes them (it fills the variables out.entropy_derivative_pressure and out.entropy_derivative_temperature), and you can then use the “latent heat” heating model to actually include the terms in the energy equation. But the latent heat material model only allows for 2 compositions, so that may not be exactly what you want either.

If you don’t want to make all of the required changes in the code, one option would also be to use the steinberger material model as the one that computes density/thermal expansivity/specific heat/compressibility when you are compositing material models. It has all of the required terms implemented, but it requires a table of material properties computed with a mineral physics software for this to work.

I hope that helps a bit, and don’t hesitate to ask if something is unclear, this is quite a lot of information and I hope I didn’t forget anything. I have only been using this feature for applications in the sublithospheric mantle, so I’ve only added the functionality to the material models designed for those sort of applications.

Best,
Juliane