diff --git a/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h b/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h
new file mode 100644
index 00000000000..4bff5d0d8da
--- /dev/null
+++ b/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h
@@ -0,0 +1,116 @@
+/*
+ Copyright (C) 2024 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file LICENSE. If not see
+ .
+*/
+
+#ifndef _aspect_material_model_rheology_compositional_plasticity_prefactors_h
+#define _aspect_material_model_rheology_compositional_plasticity_prefactors_h
+
+#include
+#include
+#include
+#include
+#include
+
+namespace aspect
+{
+ namespace MaterialModel
+ {
+ namespace Rheology
+ {
+ /**
+ * A class that handles multiplication of the plastic parameters for a given compositional
+ * field. The multiplication factors for each composition (plasticity
+ * prefactors) are also declared, parsed, and in some cases calculated in this class.
+ */
+ template
+ class CompositionalPlasticityPrefactors : public ::aspect::SimulatorAccess
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ CompositionalPlasticityPrefactors();
+
+ /**
+ * Declare the parameters this function takes through input files.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters from the parameter file.
+ */
+ void
+ parse_parameters (ParameterHandler &prm);
+
+ /**
+ * Compute the plasticity weakening factor.
+ */
+ std::pair
+ compute_weakening_factor (const MaterialModel::MaterialModelInputs &in,
+ const unsigned int volume_fraction_index,
+ const unsigned int i,
+ const double static_cohesion,
+ const double static_friction_angle,
+ const Point &position) const;
+
+ private:
+ /**
+ * The plasticity prefactors or terms used to calculate the plasticity
+ * prefactors, which are read in from the input file by the
+ * parse_parameters() function. Users can choose between different schemes.
+ * none: no change in the plasticity parameters
+ * porosity: calculate the change in the plasticity parameters due to the amount
+ * of porosity at each quadrature point.
+ * function: calculate the change in the plasticity parameters by specifying the
+ * prefactors via a function.
+ * The prefactor for a given compositional field is multiplied with a
+ * cohesions and the friction angles provided by the material model, which is then returned
+ * to the material model.
+ */
+ enum WeakeningMechanism
+ {
+ none,
+ porosity,
+ function
+ } prefactor_mechanism;
+
+ std::vector max_porosities_for_cohesion_prefactors;
+ std::vector max_porosities_for_friction_angle_prefactors;
+
+ std::vector minimum_cohesions;
+ std::vector minimum_friction_angles;
+ /**
+ * Parsed functions that specify the plasticity prefactors which must be
+ * given in the input file using the function method.
+ */
+ std::unique_ptr> prefactor_function;
+
+ /**
+ * The coordinate representation to evaluate the function for the plasticity prefactors.
+ * Possible choices are depth, cartesian and spherical.
+ */
+ Utilities::Coordinates::CoordinateSystem coordinate_system_prefactor_function;
+
+ };
+ }
+ }
+}
+#endif
diff --git a/include/aspect/material_model/rheology/visco_plastic.h b/include/aspect/material_model/rheology/visco_plastic.h
index dc9e36d05df..a46462aaedc 100644
--- a/include/aspect/material_model/rheology/visco_plastic.h
+++ b/include/aspect/material_model/rheology/visco_plastic.h
@@ -25,6 +25,7 @@
#include
#include
#include
+#include
#include
#include
#include
@@ -352,6 +353,11 @@ namespace aspect
*/
Rheology::CompositionalViscosityPrefactors compositional_viscosity_prefactors;
+ /**
+ * Object for computing the viscosity multiplied by a given prefactor term.
+ */
+ Rheology::CompositionalPlasticityPrefactors compositional_plasticity_prefactors;
+
/*
* Object for computing plastic stresses, viscosities, and additional outputs
*/
diff --git a/source/material_model/rheology/compositional_plasticity_prefactors.cc b/source/material_model/rheology/compositional_plasticity_prefactors.cc
new file mode 100644
index 00000000000..05de926b1e4
--- /dev/null
+++ b/source/material_model/rheology/compositional_plasticity_prefactors.cc
@@ -0,0 +1,228 @@
+/*
+ Copyright (C) 2024 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file LICENSE. If not see
+ .
+*/
+
+
+#include
+#include
+#include
+#include
+#include
+
+
+#include
+#include
+#include
+#include
+
+namespace aspect
+{
+ namespace MaterialModel
+ {
+ namespace Rheology
+ {
+ template
+ CompositionalPlasticityPrefactors::CompositionalPlasticityPrefactors ()
+ = default;
+
+ template
+ std::pair
+ CompositionalPlasticityPrefactors::compute_weakening_factor (const MaterialModel::MaterialModelInputs &in,
+ const unsigned int volume_fraction_index,
+ const unsigned int i,
+ const double input_cohesion,
+ const double input_friction_angle,
+ const Point &position) const
+ {
+ std::pair plasticity_parameters (input_cohesion, input_friction_angle);
+ switch (prefactor_mechanism)
+ {
+ case none:
+ {
+ return plasticity_parameters;
+ }
+ case porosity:
+ {
+ // Use the porosity to compute the cohesion and friction angle
+ // The porosity is computed in the MaterialModel::MaterialModelInputs
+ const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
+ plasticity_parameters.first *= std::max((max_porosities_for_cohesion_prefactors[volume_fraction_index] - in.composition[i][porosity_idx]) / max_porosities_for_cohesion_prefactors[volume_fraction_index], 0.0);
+ plasticity_parameters.second *= std::max((max_porosities_for_friction_angle_prefactors[volume_fraction_index] - in.composition[i][porosity_idx]) / max_porosities_for_friction_angle_prefactors[volume_fraction_index], 0.0);
+
+ plasticity_parameters.first = std::max(plasticity_parameters.first, minimum_cohesions[volume_fraction_index]);
+ plasticity_parameters.second = std::max(plasticity_parameters.second, minimum_friction_angles[volume_fraction_index] * constants::degree_to_radians);
+ return plasticity_parameters;
+ }
+ case function:
+ {
+ // Use a given function input per composition to get the prefactors
+ Utilities::NaturalCoordinate point =
+ this->get_geometry_model().cartesian_to_other_coordinates(position, coordinate_system_prefactor_function);
+
+ // we get time passed as seconds (always) but may want
+ // to reinterpret it in years
+ if (this->convert_output_to_years())
+ prefactor_function->set_time (this->get_time() / year_in_seconds);
+ else
+ prefactor_function->set_time (this->get_time());
+
+ // determine the prefactors based on position and composition
+ // The volume_fraction_index is based on the number of chemical compositional fields.
+ // However, this plugin reads a function for every compositional field, regardless of
+ // its type. Therefore we have to get the correct index.
+ // If no fields or no chemical fields are present, but only background material, the index is zero.
+ // If chemical fields are present, volume_fractions will be of size 1+n_chemical_composition_fields.
+ // The size of chemical_composition_field_indices will be one less.
+ unsigned int index = 0;
+ if (this->introspection().composition_type_exists(CompositionalFieldDescription::chemical_composition))
+ index = this->introspection().chemical_composition_field_indices()[volume_fraction_index-1];
+ double prefactor_from_function =
+ prefactor_function->value(Utilities::convert_array_to_point(point.get_coordinates()),index);
+
+ // Multiply the cohesion and the friction angle by the prefactor.
+ plasticity_parameters.first *= prefactor_from_function;
+ plasticity_parameters.second *= prefactor_from_function * constants::degree_to_radians;
+
+ return plasticity_parameters;
+ }
+ }
+ }
+
+
+ template
+ void
+ CompositionalPlasticityPrefactors::declare_parameters (ParameterHandler &prm)
+ {
+ prm.declare_entry ("Maximum porosity for cohesion prefactor", "1.0",
+ Patterns::List(Patterns::Double(std::numeric_limits::epsilon(), 1.0)),
+ "List of the maximum porosity value for calculating the prefactor for the cohesion. "
+ "entered as a volume fraction. If chosen to be 0.1 (10 percent porosity), the "
+ "cohesion will be linearly reduced from its default value when the model porosity is 0"
+ "to a maximum reduction when the model porosity is 0.1. Units: none.");
+
+ prm.declare_entry ("Maximum porosity for friction angle prefactor", "1.0",
+ Patterns::List(Patterns::Double(std::numeric_limits::epsilon(), 1.0)),
+ "List of the maximum porosity value for calculating the prefactor for the friction "
+ "angle, entered as a volume fraction. If chosen to be 0.1 (10 percent porosity), the "
+ "friction angle will be linearly reduced from its default value when the model porosity is 0"
+ "to a maximum reduction when the model porosity is 0.1. Units: none.");
+
+ prm.declare_entry ("Plasticity prefactor scheme", "none",
+ Patterns::Selection("none|function|porosity"),
+ "Select what type of plasticity multiplicative prefactor scheme to apply. "
+ "Allowed entries are 'none', 'function', and 'porosity'. 'function' allows "
+ "the cohesion and friction angle prefactors to be specified using a function. "
+ "'porosity' determines a prefactor based on the modeled porosity, this scheme "
+ "requires that the user defines a compositional field called 'porosity'. 'none' "
+ "does not modify the cohesion or the friction angle. Units: none.");
+
+ prm.declare_entry ("Minimum cohesions", "1e6",
+ Patterns::List(Patterns::Double(0.)),
+ "List of the minimum cohesions for limiting the reduction to the cohesion. Units: MPa.");
+
+ prm.declare_entry ("Minimum friction angles", "1",
+ Patterns::List(Patterns::Double(0.)),
+ "List of the minimum friction angles for limiting the reduction to the friction angle. Units: none.");
+
+ }
+
+
+
+ template
+ void
+ CompositionalPlasticityPrefactors::parse_parameters (ParameterHandler &prm)
+ {
+
+ const unsigned int n_fields = this->n_compositional_fields() + 1;
+ if (prm.get ("Plasticity prefactor scheme") == "none")
+ prefactor_mechanism = none;
+ else if (prm.get ("Plasticity prefactor scheme") == "porosity")
+ {
+ AssertThrow(this->introspection().compositional_name_exists("porosity"),
+ ExcMessage("The 'porosity' Plasticity prefactor scheme work only if "
+ "is a compositional field called porosity."));
+ prefactor_mechanism = porosity;
+
+ // Retrieve the list of chemical names
+ std::vector chemical_field_names = this->introspection().chemical_composition_field_names();
+ std::vector compositional_field_names = this->introspection().get_composition_names();
+
+ // Establish that a background field is required here
+ compositional_field_names.insert(compositional_field_names.begin(), "background");
+ chemical_field_names.insert(chemical_field_names.begin(),"background");
+
+ Utilities::MapParsing::Options options(chemical_field_names, "Maximum porosity for cohesion prefactor");
+
+ options.list_of_allowed_keys = compositional_field_names;
+ max_porosities_for_cohesion_prefactors = Utilities::MapParsing::parse_map_to_double_array (prm.get("Maximum porosity for cohesion prefactor"),
+ options);
+ max_porosities_for_friction_angle_prefactors = Utilities::MapParsing::parse_map_to_double_array (prm.get("Maximum porosity for friction angle prefactor"),
+ options);
+
+ minimum_cohesions = Utilities::MapParsing::parse_map_to_double_array (prm.get("Minimum cohesions"),
+ options);
+ minimum_friction_angles = Utilities::MapParsing::parse_map_to_double_array (prm.get("Minimum friction angles"),
+ options);
+ }
+ else if (prm.get ("Plasticity prefactor scheme") == "function")
+ {
+ prefactor_mechanism = function;
+ prm.enter_subsection("Prefactor function");
+ {
+ coordinate_system_prefactor_function = Utilities::Coordinates::string_to_coordinate_system(prm.get("Coordinate system"));
+ try
+ {
+ prefactor_function
+ = std::make_unique>(n_fields);
+ prefactor_function->parse_parameters (prm);
+ }
+ catch (...)
+ {
+ std::cerr << "FunctionParser failed to parse\n"
+ << "\t friction function\n"
+ << "with expression \n"
+ << "\t' " << prm.get("Function expression") << "'";
+ throw;
+ }
+ }
+ prm.leave_subsection();
+ }
+ else
+ AssertThrow(false, ExcMessage("Not a valid plasticity prefactor scheme"));
+ }
+ }
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace MaterialModel
+ {
+ namespace Rheology
+ {
+#define INSTANTIATE(dim) \
+ template class CompositionalPlasticityPrefactors;
+
+ ASPECT_INSTANTIATE(INSTANTIATE)
+
+#undef INSTANTIATE
+ }
+ }
+}
diff --git a/source/material_model/rheology/visco_plastic.cc b/source/material_model/rheology/visco_plastic.cc
index 94f91c241ea..106ebec2a16 100644
--- a/source/material_model/rheology/visco_plastic.cc
+++ b/source/material_model/rheology/visco_plastic.cc
@@ -351,6 +351,14 @@ namespace aspect
output_parameters.drucker_prager_parameters[j].angle_internal_friction,
in.position[i]);
+ // Step 4c: Calculate weakening factors for the cohesion and friction angles based on mechanisms
+ // other than the strain.
+ std::pair plasticity_prefactors = compositional_plasticity_prefactors.compute_weakening_factor(in, j, i,
+ output_parameters.drucker_prager_parameters[j].cohesion,
+ output_parameters.drucker_prager_parameters[j].angle_internal_friction,
+ in.position[i]);
+ output_parameters.drucker_prager_parameters[j].cohesion = plasticity_prefactors.first;
+ output_parameters.drucker_prager_parameters[j].angle_internal_friction = plasticity_prefactors.second;
// Step 5: plastic yielding
// Determine if the pressure used in Drucker Prager plasticity will be capped at 0 (default).
@@ -737,6 +745,9 @@ namespace aspect
// Variable viscosity prefactor parameters
Rheology::CompositionalViscosityPrefactors::declare_parameters(prm);
+ // Variable plasticity prefactor parameters
+ Rheology::CompositionalPlasticityPrefactors::declare_parameters(prm);
+
// Drucker Prager plasticity parameters
Rheology::DruckerPrager::declare_parameters(prm);
@@ -770,6 +781,9 @@ namespace aspect
strain_rheology.initialize_simulator (this->get_simulator());
strain_rheology.parse_parameters(prm);
+ compositional_plasticity_prefactors.initialize_simulator (this->get_simulator());
+ compositional_plasticity_prefactors.parse_parameters(prm);
+
friction_models.initialize_simulator (this->get_simulator());
friction_models.parse_parameters(prm);
diff --git a/tests/test_plastic_prefactor.prm b/tests/test_plastic_prefactor.prm
new file mode 100644
index 00000000000..c60630be9b6
--- /dev/null
+++ b/tests/test_plastic_prefactor.prm
@@ -0,0 +1,161 @@
+# A test case that checks the water fugacity viscous prefactor
+# multiplication scheme in combination with the visco plastic material model.
+# Following the flow law for wet olivine composite (diffusion and dislocation)
+# creep from Hirth & Kohlstaedt 2004 (10.1029/138GM06), the exact
+# viscosity is calculated in python through:
+# import numpy as np
+# R = 8.3144621
+# P = 1e9 Pa
+# T = 1173 K
+# edot_ii = 1e-17
+# n_disl = 3.5
+# E_disl = 520e3 J/mol
+# V_disl = 22e-6 m^3/mol
+# r_disl = 1.2
+# A_disl = 1600 / 1e6 / 1e6**n / 1e6**r 1/Pa^(4.7)/s
+#
+# n_diff = 1
+# E_diff = 375e3 J/mol
+# V_diff = 10e-6 m^3/mol
+# r_diff = 0.7
+# d = 1e-3 m
+# m = 3
+# A_diff = 2.5e7 / 1e6**n_diff / 1e6**m / 1e6**r_diff
+
+# A_H2O = 2.6e-5 1/Pa
+# activation_energy_H2O = 40e3 J/mol/K
+# activation_volume_H2O = 10e-6 m^3/mol
+#
+# M_H2O = 0.01801528 kg/mol
+# M_olivine = 0.1470027 kg/mol
+#
+# weight_fraction_h2o = 0.01
+# weight_fraction_olivine = 1 - weight_fraction_h2o
+# COH = (wt_fraction_h2o/M_H2O) / (weight_fraction_olivine/M_olivine) * 1e6
+#
+# fH2O = COH / A_H2O * np.exp((activation_energy_H2O + P*activation_volume_H2O)/(R * T))
+#
+# viscosity_disl = 0.5 * A_disl**(-1/n_disl) * edot_ii^((1-n_disl)/n_disl) * np.exp((E_disl + P*V_disl) / (n_disl*R*T)) * fH2O**(-r_disl/n_disl)
+# viscosity_diff = 0.5 * A_diff**(-1/n_diff) * d**m * np.exp((E_diff + P*V_diff) / (n_diff*R*T)) * fH2O**(-r_diff/n_diff)
+#
+# viscosity_comp = (viscosity_disl * viscosity_diff) / (viscosity_disl + viscosity_diff) = 2.62814004323392e20 Pa s
+# The viscosity should be equal to 2.6281e20 Pa s, and the material model should return this viscosity.
+
+# Global parameters
+set Dimension = 2
+set Start time = 0
+set End time = 10
+set Maximum time step = 5
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme = single Advection, iterated Stokes
+set Max nonlinear iterations = 1
+set Surface pressure = 1.e9
+
+subsection Compositional fields
+ set Names of fields = porosity
+ set Number of fields = 1
+end
+
+subsection Initial composition model
+ set Model name = function
+
+ subsection Function
+ set Function expression = if(y>=50e3, 0.05, 0)
+ end
+end
+
+# Model geometry (100x100 km, 10 km spacing)
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X repetitions = 10
+ set Y repetitions = 10
+ set X extent = 100e3
+ set Y extent = 100e3
+ end
+end
+
+# Mesh refinement specifications
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Initial global refinement = 0
+ set Time steps between mesh refinement = 0
+end
+
+# No boundary-driven deformation, but the strain-rate used by the material
+# model is set in the material model section through the parameters
+# "Reference strain rate" and "Minimum strain rate".
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = bottom, top, left, right
+end
+
+# Initial temperature conditions
+subsection Initial temperature model
+ set Model name = function
+
+ subsection Function
+ set Function expression = 1173
+ end
+end
+
+subsection Material model
+ set Model name = visco plastic
+
+ subsection Visco Plastic
+ # Composite creep, and choose the Hirth & Kohlstaedt 2004 olivine hydration
+ # viscosity prefactor scheme.
+ set Viscous flow law = composite
+ set Reference strain rate = 1e-17
+ set Minimum strain rate = 1e-17
+
+ # Dislocation creep parameters
+ # The water fugacity exponent is 1.2, but we need to divide by
+ # the stress exponent, so we input r/n=0.34285714285714286.
+ set Prefactors for dislocation creep = 1.0095317511683098e-25
+ set Stress exponents for dislocation creep = 3.5
+ set Activation energies for dislocation creep = 520e3
+ set Activation volumes for dislocation creep = 22e-6
+ set Water fugacity exponents for dislocation creep = 0.34285714285714286
+
+ # Diffusion creep parameters
+ # The same approach is required as for the dislocation parameters,
+ # but the stress exponent for diffusion creep is 1.
+ set Prefactors for diffusion creep = 1.5773933612004842e-21
+ set Stress exponents for diffusion creep = 1.0
+ set Activation energies for diffusion creep = 375e3
+ set Activation volumes for diffusion creep = 10e-6
+ set Grain size = 1e-3
+ set Water fugacity exponents for diffusion creep = 0.7
+
+ set Cohesions = 100e6
+ set Angles of internal friction = 30
+ set Plasticity prefactor scheme = porosity
+ set Maximum porosity for friction angle prefactor = 0.
+ set Maximum porosity for cohesion prefactor = 0.1
+
+ set Minimum cohesions = 10e6
+ set Minimum friction angles = 25
+ end
+end
+
+# Gravity model
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 0.0
+ end
+end
+
+# Post processing
+subsection Postprocess
+ set List of postprocessors = material statistics, visualization
+
+ subsection Visualization
+ set Interpolate output = false
+ set Time between graphical output = 0
+ set List of output variables = material properties, named additional outputs
+ set Output format = vtu
+ end
+end