Melt Production/ Migration

ASPECT
1

Hello!
I’m a grad student exploring/ modeling melt migration in an early continental rift setting. I’m currently setting up a buoyancy driven viscous flow to model rifting initiation and partial melting using ASPECT. To my understanding, melt production is calculated using the Katz et al., 2003, and Sobolev et al., 2011 for peridotite and pyroxenite, respetively, and then Keller et al., 2013 formulations, for viscoelastic conditions is that correct?
I’m exploring my options in parameterizing and was curious if there a way to modify melt parameters so that they don’t exactly follow the Katz et al., 2003 formulation i.e. adding an opx phase? Also are we limited to tracking melt through a porosity composition field or is there a way to use particles to track melt.

Cordially,
Carlos

2

Hi Carlos,

Juliane Dannberg will likely weigh in regarding the different melting models, so I will not touch on that.

Regarding the viscoelastic (+ plastic) implementation, we are working on the Keller et al. formulation, but do not have a version of that ready for general use yet and it is likely at least 6-12 months out.

However, it would be possible to have a VEP rheology with melt, where the fluid pressure does not influence the yielding condition or viscoelastic formulation.

If you would like to go down the latter path, we should discuss this in detail with others at one of the bi-weekly user/developer meetings.

Cheers,
John

3

Hi Carlos,

great to hear that you want to use ASPECT for modeling melting and melt migration!
Because this is such a complex problem, it’s not that simple to answer your question. So let me address your questions separately:

(1) Melting model.
The melting model depends on the material model that you use. There is one that uses the Katz et al., 2003, and Sobolev et al., 2011 melting models (this is the latent heat melt material model), but this model does not include melt migration. There is also the melt simple model, which has both melting and melt migration and uses the (anhydrous part of the) Katz 2003 model. However, there is a problem with using that parametrizaton for models with melt migration: it does not predict the melt composition. Without that, as soon as melt migrates away from where it was generated, it is not clear what the bulk composition is, and so the model does not predict where melt should be freezing again. The ASPECT material model implements a workaround for this, but it’s not thermodynamically consistent.

This is a limitation of the original melting model (not the implementation in ASPECT), and commonly, people who use models of melt migration use other melting models, which are simpler in a sense that they don’t track the different mineral phases in the solid (the the ‘kink’ in the Katz 2003 model wouldn’t be included).

You can of course modify the melting parametrization in any way you want (like adding an opx phase), but these are some of the limitations to keep in mind.

(2) The mechanical model for melt migration.
ASPECT has the option to solve the McKenzie equations for two-phase flow. At the moment, this is implemented in material models with a viscous rheology. It is not too difficult to modify other material models to include melt migration, but as soon as they have a plastic or elastic rheology component, additional terms are required in the equations, and it also becomes numerically really hard to solve the equations (for example due to the hydrofracturing terms related to the interaction of melt and brittle failure). Several members of the ASPECT team are working on different parts of these problems, but this is just a really complex problem. So at the moment, it is not possible to use both melt migration and visco-elasto-plastic deformation in the same model (but let us know if you want to help with the implementation!).

(3) Melt and particles.
It is possible to move particles with the melt velocity, but the porosity equation would still be solved using a compositional field. This is because it has extra terms in the equations that are different from a normal compositional field (but means you can track the melt composition on particles). But the other problem with using particles here is that you may not want to have the melt particles everywhere in the model, but only where melt is present, so you would have to add a particle generator routine that adds/removes particles based on melting or freezing, and it is not straightforward how this should look like.

So my question to you would be: What is the purpose of your model, and what specific components/physical processes do you need for that? For example, is it important to see where the melt migrates to, or only where it is being generated?

Best,
Juliane

1 Like
4

Hi Carlos, John and Juliane,

Thanks for the interesting discussion, and I’d love to learn more about all of these topics. John also mentions the bi-weekly user/developer meeting. I’d be interested in attending such meetings (perhaps not every one of them, but whenever I would have time). If these are open to all and suitable for those with not much too developer experience (but eager to learn), could you please give some information about those?

Thanks,

Jeroen

5

Hi Jeroen,

Absolutely, the user meetings are open to anyone interested in participating! Participation includes just listening to see what others are doing and keeping in touch, and there is no requirement to present anything.

The announcements for the bi-weekly meetings are done on this forum, with updates send the Tuesday before each meeting:

The link above also contains a zoom link and a link to a google doc, where anyone can add discussion topics.

It would be great to have more users attend!

Cheers,
John

6

Hi,

I am also interested in the modelling of melt migration using ASPECT.
I tried to set up a simple box model with the “melt simple” material heated from below. If I prescribe a fixed temperature at the top and bottom boundaries, the results look reasonable. However, it would be much more convenient for me to prescribe a (non-zero) heat flux at the bottom boundary. I tried the “Boundary heat flux model”, but it apparently doesn’t work when melt migration is modelled.

Is it a problem to prescribe the heat flux BC in a model with melt migration? Or it just hasn’t been implemented (yet)?

Best wishes,
Petra

7

Hi Petra,

it’s great that you’re interested in using melt migration in ASPECT!
can you be a bit more precise what you mean by “it apparently doesn’t work”? Do you get an error message, or do the results look like there is no prescribed heat flux?

In principle, there shouldn’t be a problem with using both features together, they just both add/replace terms in the energy equation and I haven’t tested if both features work together.

Best regards,
Juliane

8

Hi,
thanks for your quick reply!
The model runs without error messages, but there is no heating at the bottom boundary, even though the model statistics and output in Paraview (heat flux map) show the correct value of the heat flux. The change happens when I switch on the melt transport:
set Include melt transport = true
In the same model setup, when I switch it off, the heating at the bottom boundary is clearly visible.
Petra

9

Hi Petra,

this sounds like a bug in ASPECT, thank you for letting us know! I’ll have a look to see if we can fix it. I talked to @gassmoeller and he already had an idea what might be the problem. We’ll let you know how this goes.

Best regards,
Juliane

10

OK, thanks a lot Juliane!
Petra

11

This should be fixed now (see this pull request for more details). Can you update your aspect to the current developer version and see if that works for your model?

Best regards,
Juliane

12

Hi Juliane, hi Rene,
I’ve tested the new version and it works perfectly!
Now it is time for me to play with the models;)
Thank you very much!
With best regards,
Petra

13

Hi Juliane & John,

Thanks for your thorough answers! I’ve been working on my melt model and wanted to follow up with a few questions. I wanted to know how to best go about discretizing a high-density anomaly in the lower lithosphere to initiate a negatively-buoyant lithospheric drip. The purpose of the model would be to analyze where melt would form as this falling drip meets thermally buoyant mantle material. My main problem right now is in adding another compositional field, how do I prescribe a density for a field in the ‘melt simple’ model?

Cordially,
Carlos G

14

Hi Carlos,

the ‘melt simple’ model currently only takes into account density changes by one compositional field, the one that tracks the depletion (with the name “peridotite”). So in principle, you could use just this one compositional field called peridotite, and use the composition initial conditions to set the composition to a negative value where you want your dense drip to be. You would also have to make sure you have set the “Depletion density change” to a negative value (indicating lower densities for depleted material). I haven’t actually tested this, so I don’t know how the melt parametrization would handle this (the Katz 2003 parametrization is for batch melting of peridotite, so it does not cover this type of composition).

But the other question I would think about is this: Do you only want to see where melt would form? In this case you can use any material model, and simply use the “melt fraction” postprocessor.

The melt material models are only useful if you want to know how melt migrates. I’ve actually run a model that self-consistently develops lithospheric drips (see for example here) using the melt global material model as a proof-of-concept for this paper, so that material model is actually designed for models just like the one you describe. However, this has the disadvantage that the melting parametrization is really simplified.

Best,
Juliane

15

Hi Carlos,

I agree with Juliane’s suggestions above if you can address the questions in hand with only a mantle rheology or you just want to look at where melt forms (i.e., not melt migration).

If you want to include both a complicated lithospheric structure + visco-plastic or viscoelastic-plastic rheology and melt migration, you will need to do some additional coding. If this is the case, a new material model plugin that derives from an existing material model is likely the best path forward and we can discuss this in more detail.

Cheers,
John

16

Hi John,

I am trying to do what you are suggesting, including both visco-plastic rheology with melt migration. To that end, I tried to merge the code from one of the tests called rising_melt_blob.cc with visco-plastic rheology code in a new plugin. I added the following part from there

 template <int dim>
    class MeltingRate:
    public MaterialModel::MeltInterface<dim>, public ::aspect::SimulatorAccess<dim>
  {
      virtual bool is_compressible () const
      {
        return false;
      }

      virtual double reference_viscosity () const
      {
        return 5e20;
      }

      virtual double reference_darcy_coefficient () const
      {
        return 1e-8 * std::pow(0.01, 3.0) / 10.0;
      }

      virtual void evaluate(const typename MaterialModel::Interface<dim>::MaterialModelInputs &in,
                            typename MaterialModel::Interface<dim>::MaterialModelOutputs &out) const
      {

        // fill melt outputs if they exist
        aspect::MaterialModel::MeltOutputs<dim> *melt_out = out.template get_additional_output<aspect::MaterialModel::MeltOutputs<dim> >();

        if (melt_out != nullptr)
          {
            const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
            for (unsigned int i=0; i<in.n_evaluation_points(); ++i)
              {
                const double porosity = in.composition[i][porosity_idx];

                melt_out->compaction_viscosities[i] = 5e20 * 0.05/std::max(porosity,0.00025);
                melt_out->fluid_viscosities[i]= 10.0;
                melt_out->permeabilities[i]= 1e-8 * std::pow(porosity,3) * std::pow(1.0-porosity,2);
                melt_out->fluid_densities[i]= 2000.0;
                melt_out->fluid_density_gradients[i] = Tensor<1,dim>();
              }
          }
      }

  };

I declared an additional composition called porosity and tracked it with melt particle as is done in particle_melt_advection.prm from tests.

While compiling this new plugin I got no errors but running the model gives this error


----------------------------------------------------
Exception 'ExcMessage("You have requested to convert a plugin of type <" + boost::core::demangle(typeid(PluginType).name()) + "> into type <" + boost::core::demangle(typeid(TestType).name()) + ">, but this cast cannot be performed.")' on rank 0 on processing: 

--------------------------------------------------------
An error occurred in line <87> of file </home/mainakchakraborty/bin/aspect/include/aspect/plugins.h> in function
    TestType& aspect::Plugins::get_plugin_as_type(PluginType&) [with TestType = const aspect::MaterialModel::MeltInterface<2>; PluginType = const aspect::MaterialModel::Interface<2>]
The violated condition was: 
    plugin_type_matches<TestType>(object)
Additional information: 
    You have requested to convert a plugin of type
    <aspect::MaterialModel::Interface<2>> into type
    <aspect::MaterialModel::MeltInterface<2>>, but this cast cannot be
    performed.
--------------------------------------------------------

Aborting!
----------------------------------------------------

Is there another way we can derive from the existing melt migration code. Please let me know.

Thanks,
Mainak