# Plastic strain in reactive fluid transport rheology

**URL:** <https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996>\
**Category:** ASPECT\
**Created:** [March 11, 2025, 11:29am UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996 "2025-03-11T11:29:07Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![fgueyd](https://avatars.discourse-cdn.com/v4/letter/f/ecc23a/32.png) [@fgueyd](https://community.geodynamics.org/u/fgueyd)\
**Post date:** [March 11, 2025, 11:29am UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/1 "2025-03-11T11:29:07Z")

</div>

Dear all,

I used the material model “reactive fluid transport “ based on the cookbook example developed by Daniel Douglas ([Simplified Subduction Model with Parameterized Solid-Fluid Reactions — ASPECT 3.1.0-pre](https://aspect-documentation.readthedocs.io/en/latest/user/cookbooks/cookbooks/tian_parameterization_kinematic_slab/doc/tian_parameterization_kinematic_slab.html)).

The base model is visco-plastic

I wish to track the plastic strain through time, as it is done in the classical visco-plastic rheology

I define a compositional field “plastic\_strain” and define a Strain weakening mechanism so that the plastic strain is computed

The model run perfectly but the plastic strain is zero

It seems that when “reactive fluid transport” material model is used, even if the base model is visco-plastic, the plastic strain is not computed. I saw this message in the log  
“Skipping plastic\_strain composition solve because RHS is zero.”

Could you please let me know how to fix this problem

I attached hereafter two prm files  
-one with a simple extension test with only visco-plastic rheology (plastic strain is well computed)  
-one with the same parameters except the use of “reactive fluid transport” rheology (plastic strain is not computed)

Thank you for your feedbacks

Frédéric

* * *

## -prm with a simple extension test with only visco-plastic rheology (plastic strain is well computed)

set Dimension = 2  
set Use years in output instead of seconds = true  
set Output directory = strain-VP-iterated  
set End time = 0.5e6 #years

set Nonlinear solver scheme = iterated Advection and Stokes  
set Nonlinear solver tolerance = 1e-1 # normaly 1e-5  
set Max nonlinear iterations = 5 # 100 normaly

set Adiabatic surface temperature = 1600  
#set Use operator splitting = true

subsection Geometry model  
set Model name = box  
subsection Box  
set X repetitions = 4  
set Y repetitions = 2  
set X extent = 100e3  
set Y extent = 50e3  
end  
end

subsection Mesh deformation  
set Mesh deformation boundary indicators = top: free surface  
subsection Free surface  
set Surface velocity projection = vertical  
end  
set Additional tangential mesh velocity boundary indicators = left, right  
end

set Pressure normalization = no

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

subsection Initial temperature model  
set Model name = function  
subsection Function  
set Variable names = x,y,t  
set Function constants = Ttop=273, Tbot=1473, ysize=50e3  
set Function expression = Ttop+(Tbot-Ttop)\*(ysize-y)/ysize  
end  
end

subsection Boundary temperature model  
set Fixed temperature boundary indicators = top, bottom  
set List of model names = function  
subsection Function  
set Function constants = Ttop=273, Tbot=1473, ysize=50e3  
set Function expression = (y\>ysize/2 ? Ttop : Tbot)  
end  
end

subsection Boundary velocity model  
set Prescribed velocity boundary indicators = bottom y: function, left x: function, right x: function  
subsection Function  
set Variable names = x,y  
set Function constants = v=0.06, xsize=100e3, ysize=50e3  
set Function expression = (x\<xsize/2 ? -v : v); (y\<ysize/2 ? ysize/xsize_v : -ysize/xsize_v)  
end  
end

subsection Compositional fields  
set Number of fields = 4  
set Names of fields = peridotite, porosity, plastic\_strain, seed  
set Compositional field methods = field, field, field, field  
end

subsection Initial composition model  
set List of model names = function  
subsection Function  
set Variable names = x,y  
set Function constants = xsize=100e3, ypos=30e3, rad=5e3, x1=40e3,x2=60e3,y1=5e3  
set Function expression = 0; if(x\>x1 && x\<x2 && y\<y1, 0.00, 0);0; if((x-xsize/2)^2+(y-ypos)^2\<rad^2, 1, 0)  
end  
end

subsection Boundary composition model  
set Fixed composition boundary indicators = top, bottom, left, right  
set List of model names = initial composition  
end

subsection Heating model  
set List of model names = compositional heating

subsection Compositional heating  
set Use compositional field for heat production averaging = 0  
set Compositional heating values = 0.0  
end  
end

subsection Material model  
set Model name = visco plastic  
set Material averaging = harmonic average only viscosity

subsection Visco Plastic  
set Densities = 3300 # kg/m3  
set Thermal diffusivities = 1e-6  
set Heat capacities = 1000  
set Thermal expansivities = 0

```
set Viscous flow law = dislocation

set Prefactors for dislocation creep = 2.7e-18, 1,1,1, 1e-15
set Activation energies for dislocation creep = 480e3,0,0,0,480e3
set Activation volumes for dislocation creep = 1.1e-5,0,0,0,1.1e-5
set Stress exponents for dislocation creep = 3.5,1,1,1,3.5

set Angles of internal friction = 10.0, 10.,10, 1, 9.0 # degrees
set Cohesions = 10e6 # Pa

set Minimum viscosity = 1e19 # Pa s
set Maximum viscosity = 1e25 # Pa s
set Strain weakening mechanism = plastic weakening with plastic strain only
set Start plasticity strain weakening intervals = 0.1
set End plasticity strain weakening intervals = 1
set Cohesion strain weakening factors = 0.25
set Friction strain weakening factors = 0.25

```

end  
end

subsection Mesh refinement  
set Initial global refinement = 2 # normally 4  
set Initial adaptive refinement = 2  
set Time steps between mesh refinement = 1  
set Strategy = composition, strain rate  
set Refinement criteria merge operation = max  
set Refinement criteria scaling factors = 1, 2  
end

subsection Melt settings  
#set Include melt transport = true  
end

subsection Postprocess  
set List of postprocessors = visualization  
subsection Visualization  
set Time between graphical output = 1  
set Output mesh velocity = true  
set List of output variables = material properties, strain rate, heating  
subsection Material properties  
set List of material properties = density, viscosity  
end  
end  
end

* * *

## -prm with the same parameters except the use of “reactive fluid transport” rheology (plastic strain is not computed)

set Dimension = 2  
set Use years in output instead of seconds = true  
set Output directory = strain-ReactiveFluidTransport-iterated  
set End time = 0.5e6 #years

set Nonlinear solver scheme = iterated Advection and Stokes  
set Nonlinear solver tolerance = 1e-1 # normaly 1e-5  
set Max nonlinear iterations = 5 # 100 normaly

set Adiabatic surface temperature = 1600  
#set Use operator splitting = true

subsection Geometry model  
set Model name = box  
subsection Box  
set X repetitions = 4  
set Y repetitions = 2  
set X extent = 100e3  
set Y extent = 50e3  
end  
end

subsection Mesh deformation  
set Mesh deformation boundary indicators = top: free surface  
subsection Free surface  
set Surface velocity projection = vertical  
end  
set Additional tangential mesh velocity boundary indicators = left, right  
end

set Pressure normalization = no

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

subsection Initial temperature model  
set Model name = function  
subsection Function  
set Variable names = x,y,t  
set Function constants = Ttop=273, Tbot=1473, ysize=50e3  
set Function expression = Ttop+(Tbot-Ttop)\*(ysize-y)/ysize  
end  
end

subsection Boundary temperature model  
set Fixed temperature boundary indicators = top, bottom  
set List of model names = function  
subsection Function  
set Function constants = Ttop=273, Tbot=1473, ysize=50e3  
set Function expression = (y\>ysize/2 ? Ttop : Tbot)  
end  
end

subsection Boundary velocity model  
set Prescribed velocity boundary indicators = bottom y: function, left x: function, right x: function  
subsection Function  
set Variable names = x,y  
set Function constants = v=0.06, xsize=100e3, ysize=50e3  
set Function expression = (x\<xsize/2 ? -v : v); (y\<ysize/2 ? ysize/xsize_v : -ysize/xsize_v)  
end  
end

subsection Compositional fields  
set Number of fields = 4  
set Names of fields = peridotite, porosity, plastic\_strain, seed  
set Compositional field methods = field, field, field, field  
end

subsection Initial composition model  
set List of model names = function  
subsection Function  
set Variable names = x,y  
set Function constants = xsize=100e3, ypos=30e3, rad=5e3, x1=40e3,x2=60e3,y1=5e3  
set Function expression = 0; if(x\>x1 && x\<x2 && y\<y1, 0.00, 0);0; if((x-xsize/2)^2+(y-ypos)^2\<rad^2, 1, 0)  
end  
end

subsection Boundary composition model  
set Fixed composition boundary indicators = top, bottom, left, right  
set List of model names = initial composition  
end

subsection Heating model  
set List of model names = compositional heating

subsection Compositional heating  
set Use compositional field for heat production averaging = 0  
set Compositional heating values = 0.0  
end  
end

subsection Material model  
set Model name = reactive fluid transport  
set Material averaging = harmonic average  
subsection Reactive Fluid Transport Model  
set Base model = visco plastic  
set Reference fluid density = 3200 # density of water  
set Shear to bulk viscosity ratio = 0.1  
set Reference fluid viscosity = 1.  
set Reference permeability = 1e-6  
set Exponential fluid weakening factor = 5  
set Fluid compressibility = 0  
set Fluid reaction time scale for operator splitting = 5e4  
set Fluid-solid reaction scheme = katz2003  
subsection Katz 2003 model  
set A1 = 941.32  
set A2 = 1.3463e-7  
set A3 = -5.3158e-18  
end  
end  
subsection Visco Plastic  
set Densities = 3300 # kg/m3  
set Thermal diffusivities = 1e-6  
set Heat capacities = 1000  
set Thermal expansivities = 0

```
set Viscous flow law = dislocation

set Prefactors for dislocation creep = 2.7e-18, 1,1,1, 1e-15
set Activation energies for dislocation creep = 480e3
set Activation volumes for dislocation creep = 1.1e-5
set Stress exponents for dislocation creep = 3.5

set Angles of internal friction = 10.0, 10.,10, 1, 9.0 # degrees
set Cohesions = 10e6 # Pa

set Minimum viscosity = 1e19 # Pa s
set Maximum viscosity = 1e25 # Pa s
set Strain weakening mechanism = plastic weakening with plastic strain only
set Start plasticity strain weakening intervals = 0.1
set End plasticity strain weakening intervals = 1
set Cohesion strain weakening factors = 0.25
set Friction strain weakening factors = 0.25

```

end  
end

subsection Mesh refinement  
set Initial global refinement = 2 # normally 4  
set Initial adaptive refinement = 2  
set Time steps between mesh refinement = 1  
set Strategy = composition, strain rate  
set Refinement criteria merge operation = max  
set Refinement criteria scaling factors = 1, 2  
end

subsection Melt settings  
#set Include melt transport = true  
end

subsection Postprocess  
set List of postprocessors = visualization  
subsection Visualization  
set Time between graphical output = 1  
set Output mesh velocity = true  
set List of output variables = material properties, strain rate, heating, melt fraction #, melt material properties  
subsection Material properties  
set List of material properties = density, viscosity  
end  
end  
end

---

<div class="post-metadata">

**Author:** ![daniel.douglas](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/daniel.douglas/32/2126_2.png) [@daniel.douglas](https://community.geodynamics.org/u/daniel.douglas)\
**Post date:** [March 11, 2025, 8:49pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/2 "2025-03-11T20:49:46Z")

</div>

Dear Frédéric

Thanks for posting on the forum and happy to see that this cookbook example might be of use to some people!! I haven’t personally tried running a model with the `Reactive Fluid Transport Model` while tracking `plastic_strain`, but I think that this should work. From looking at your .prm file, I suspect that the reason it is saying that the plastic\_strain is 0 is because you are beginning with a `plastic_strain` that is equal to 0 everywhere:

```auto
subsection Compositional fields
  set Number of fields = 4
  set Names of fields = peridotite, porosity, plastic_strain, seed
  set Compositional field methods = field, field, field, field
end

subsection Initial composition model
  set List of model names = function
  subsection Function
    set Variable names = x,y
    set Function constants = xsize=100e3, ypos=30e3, rad=5e3, x1=40e3,x2=60e3,y1=5e3
    set Function expression = 0; if(x>x1 && x<x2 && y<y1, 0.00, 0);0; if((x-xsize/2)^2+(y-ypos)^2<rad^2, 1, 0)
  end
end

```

If defining an initial `plastic_strain` does not help to overcome the issue, I would try outputting `named additional outputs`:

```auto
subsection Visualization
  set Time between graphical output = 1
  set Output mesh velocity = true
  set List of output variables = material properties, strain rate, heating, named additional outputs

```

This will tell you where in the model that there is actually plastic yielding occurring, perhaps the model is not actually yielding anywhere?

I hope this helps!  
Cheers,  
Daniel

---

<div class="post-metadata">

**Author:** ![fgueyd](https://avatars.discourse-cdn.com/v4/letter/f/ecc23a/32.png) [@fgueyd](https://community.geodynamics.org/u/fgueyd)\
**Post date:** [March 11, 2025, 9:22pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/3 "2025-03-11T21:22:30Z")

</div>

Hi Daniel  
thank you for your prompt answer

I tried these two solutions but it is not working  
I think that simply the plastic strain is not calculated  
If I use ViscoPlastic instead of Reactive Fluid Transport (for materail model), plastic strain is calculated.

So I guess that in Reactive Fluid Transport rheology, the call to strain dependant module (or something like that) is not done

Any guess?

Frédéric

---

<div class="post-metadata">

**Author:** ![fgueyd](https://avatars.discourse-cdn.com/v4/letter/f/ecc23a/32.png) [@fgueyd](https://community.geodynamics.org/u/fgueyd)\
**Post date:** [March 19, 2025, 4:32pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/4 "2025-03-19T16:32:39Z")

</div>

Dear all and in particular aspect developers/experts

Do you have any ideas if it is possible to track plastic strain (or total strain) in reactive fluid transport Material Model?

Many thanks for your help

Best

Frédéric

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [March 19, 2025, 5:07pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/5 "2025-03-19T17:07:49Z")

</div>

@fgueyd - Sorry for not chiming into this discussion earlier.

I’m not entirely sure what the issue is with your setup, but I did run a quick test where viscous strain weakening was added an existing test (`reactive_fluid_transport_zero_solubility.prm`) that uses the reactive fluid transport material model with the visco plastic material model as the base model.

Viscous strain does accumulate in this model after 1 time step.

A few questions and requests:

1. Using the attached PRM as a template, can you try replacing plastic strain with viscous strain in your model and see if it accumulates?
2. Can you confirm that plastic yielding is occurring in your model?
3. This should not be an issue due to internal logic within the code, but it might be worth specifying the [field type](https://aspect-documentation.readthedocs.io/en/latest/parameters/Compositional_20fields.html#parameter-name-types-of-fields) for each composition you are tracking (see also attached PRM).
4. Develop a small test case, perhaps building on the attached PRM, where yielding should occur and you know exactly how much plastic strain should develop in each time step? This can be done by setting` Minimum strain rate` equal to the `Reference strain rate`.

Cheers,  
John  
[reactive\_fluid\_transport\_zero\_solubility\_viscous\_strain.prm](https://community.geodynamics.org/uploads/short-url/9wPX5o5Y3dranT0PrqWtnhu622T.prm) (4.7 KB)

---

<div class="post-metadata">

**Author:** ![jdannberg](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/jdannberg/32/26_2.png) [@jdannberg](https://community.geodynamics.org/u/jdannberg)\
**Post date:** [March 19, 2025, 8:15pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/6 "2025-03-19T20:15:15Z")

</div>

Hi everyone,

I suspect the problem is the following: Frédéric, you use the katz2003 fluid reaction scheme, and checking what that scheme does, it sets all reaction\_terms to zero (that’s what happens when `katz2003_model.calculate_reaction_rate_outputs(in, out)` is called in line [216](https://github.com/geodynamics/aspect/blob/5f73a331118c291b55f4ec32a6ef702f8aca1995/source/material_model/reactive_fluid_transport.cc#L216C11-L216C67)). The corresponding line in katz2003\_mantle\_melting.cc is [here](https://github.com/geodynamics/aspect/blob/5f73a331118c291b55f4ec32a6ef702f8aca1995/source/material_model/reaction_model/katz2003_mantle_melting.cc#L281). However, these reaction\_terms also determine how the finite strain changes in the strain dependent rheology. But since the evaluate() function of the base model is called first (this is where the viscoplastic model fills the reaction\_terms) and then katz2003\_model.calculate\_reaction\_rate\_outputs is called afterwards, it would overwrite these terms, setting them to zero instead. This would explain why you do not see any change in the plastic strain.

So I suspect you could fix your problem by modifying katz2003\_mantle\_melting.cc so that it no longer sets all reaction\_terms to zero.

Best,  
Juliane

---

<div class="post-metadata">

**Author:** ![fgueyd](https://avatars.discourse-cdn.com/v4/letter/f/ecc23a/32.png) [@fgueyd](https://community.geodynamics.org/u/fgueyd)\
**Post date:** [March 20, 2025, 5:44pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/7 "2025-03-20T17:44:57Z")

</div>

dear all  
thank your very much for your inputs. Now it is working  
@jbnaliboff : your test indeed worked but not with katz2003  
@jdannberg : perfect, it is working perfectly. I was looking at this .cc file but I am not familiar with c++

Best  
Frédéric

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [March 21, 2025, 9:35pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/8 "2025-03-21T21:35:04Z")

</div>

@fgueyd - I opened a [pull request](https://github.com/geodynamics/aspect/pull/6265/files) with the suggested fix from @jdannberg.

A quick test (see attached) confirmed plastic strain now accumulates, but it would be great to see if the branch works for your model (would you mind testing this?).

Cheers,  
John

[viscoplastic\_melt\_fraction\_plastic\_strain.prm](https://community.geodynamics.org/uploads/short-url/zo4KT4FAMZ5605BR4ocZla6hmv1.prm) (4.4 KB)

---

<div class="post-metadata">

**Author:** ![fgueyd](https://avatars.discourse-cdn.com/v4/letter/f/ecc23a/32.png) [@fgueyd](https://community.geodynamics.org/u/fgueyd)\
**Post date:** [March 22, 2025, 8:15am UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/9 "2025-03-22T08:15:39Z")

</div>

@jbnaliboff I am testing it. It is working but it seems that the melting (e.g porosity) is occuring in the same way as previously (e.g. without the modification proposed by @jdannberg). I am conducing further tests to understand this. I will let you know  
Cheers  
Frédéric

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [March 22, 2025, 3:46pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/10 "2025-03-22T15:46:46Z")

</div>

> @fgueyd - am testing it. It is working but it seems that the melting (e.g porosity) is occuring in the same way as previously (e.g. without the modification proposed by @jdannberg).

Yes, I found one other place that needed a similar conditional check, can you check the updated PR now? Note that these change introduced some additional complications to consider (see discussion in PR), and this may take a bit of time to merge and fully test. Please proceed with caution when using these changes until the PR is fully merged and tested.

---

<div class="post-metadata">

**Author:** ![fgueyd](https://avatars.discourse-cdn.com/v4/letter/f/ecc23a/32.png) [@fgueyd](https://community.geodynamics.org/u/fgueyd)\
**Post date:** [March 22, 2025, 4:49pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/11 "2025-03-22T16:49:05Z")

</div>

thank you John. I will check the PR and test it

---

<div class="post-metadata">

**Author:** ![gassmoeller](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/gassmoeller/32/11_2.png) [@gassmoeller](https://community.geodynamics.org/u/gassmoeller)\
**Post date:** [March 27, 2025, 3:46pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/12 "2025-03-27T15:46:07Z")

</div>

Hi @fgueyd,  
I just merged a potential solution by @jdannberg into the ASPECT main branch ([fill reaction terms in material model, not reaction model by jdannberg · Pull Request #6266 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/pull/6266)). It is an improvement of the code style either way, but could you check if your problem is fixed in the latest development version and report back to us? Hopefully the problem should be fixed now.

Thanks,  
Rene

---

<div class="post-metadata">

**Author:** ![fgueyd](https://avatars.discourse-cdn.com/v4/letter/f/ecc23a/32.png) [@fgueyd](https://community.geodynamics.org/u/fgueyd)\
**Post date:** [April 1, 2025, 6:41pm UTC](https://community.geodynamics.org/t/plastic-strain-in-reactive-fluid-transport-rheology/3996/13 "2025-04-01T18:41:29Z")

</div>

Hi @gassmoeller, I will let you know as soon as possible (lots of courses this week..)