Skip to content

[WIP] Tectonic-dependent dynamic friction weakening #6541

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
+132 −22
Open

[WIP] Tectonic-dependent dynamic friction weakening #6541

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a14991d
add differential strain rate weakening
Minerallo Jun 18, 2025
70ae7e4
add Differential strain rate weakening case
Minerallo Jun 19, 2025
40cbb75
Corrected variable type
Minerallo Jun 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 28 additions & 3 deletions include/aspect/material_model/rheology/friction_models.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@

You should have received a copy of the GNU General Public License
along with ASPECT; see the file LICENSE. If not see
<http://www.gnu.org/licenses/>.
<http://www.gnu.org/licenses/>.
*/

#ifndef _aspect_material_model_rheology_friction_models_h
Expand Down Expand Up @@ -43,11 +43,15 @@ namespace aspect
*
* For the type 'dynamic friction', the friction angle is rate dependent following
* Equation 13 from \\cite{van_dinther_seismic_2013}.
*
* For the type 'differential dynamic friction', the dynamic angle depends on whether
* local flow is convergent or divergent.
*/
enum FrictionMechanism
{
static_friction,
dynamic_friction,
differential_dynamic_friction,
function
};

Expand Down Expand Up @@ -77,7 +81,8 @@ namespace aspect
compute_friction_angle(const double current_edot_ii,
const unsigned int volume_fraction_index,
const double static_friction_angle,
const Point<dim> &position) const;
const Point<dim> &position,
const SymmetricTensor<2, dim> &strain_rate) const;

/**
* A function that returns the selected type of friction dependence.
Expand All @@ -88,7 +93,7 @@ namespace aspect
private:
/**
* Select the mechanism to be used for the friction dependence.
* Possible options: static friction | dynamic friction | function
* Possible options: static friction | dynamic friction | differential_dynamic_friction | function |
*/
FrictionMechanism friction_mechanism;

Expand Down Expand Up @@ -127,6 +132,26 @@ namespace aspect
* Possible choices are depth, cartesian and spherical.
*/
Utilities::Coordinates::CoordinateSystem coordinate_system_friction_function;

/**
* Dynamic angles of internal friction that are used at high strain rates in converging regions,
* when using differential dynamic friction.
*/
std::vector<double> dynamic_angles_of_internal_friction_for_convergence;
/**
* Dynamic angles of internal friction that are used at high strain rates in diverging regions,
* when using differential dynamic friction.
*/
std::vector<double> dynamic_angles_of_internal_friction_for_divergence;

/**
* Thresholds on the strain rate trace to classify convergent regimes.
*/
double convergence_threshold;
/**
* Thresholds on the strain rate trace to classify divergent regimes.
*/
double divergence_threshold;
};
}
}
Expand Down
120 changes: 102 additions & 18 deletions source/material_model/rheology/friction_models.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ namespace aspect
compute_friction_angle(const double current_edot_ii,
const unsigned int volume_fraction_index,
const double static_friction_angle,
const Point<dim> &position) const
const Point<dim> &position,
const SymmetricTensor<2, dim> &strain_rate) const
{

switch (friction_mechanism)
Expand Down Expand Up @@ -76,6 +77,46 @@ namespace aspect
"The friction angle should be smaller than 1.6 rad."));
return dynamic_friction_angle;
}
case differential_dynamic_friction:
{
// Calculate effective steady-state friction coefficient.
// This variant of the dynamic friction model allows the minimum friction angle
// to differ between convergent and divergent settings. It reflects the idea that
// different tectonic environments exhibit different weakening behavior — for example,
// lower friction in subduction zones due to fluids, or in ridges due to melt.
// The dynamic angle used at high strain rates is chosen based on the local flow
// (converging or diverging), while the static angle applies at low strain rates.
//
// The transition is smooth, following:
// mu = mu_d + (mu_s - mu_d) / (1 + (strain_rate / strain_rate_c)^X)
//
// where mu_d is selected based on convergence/divergence.
// A 3-zone dynamic friction model where the minimum friction angle varies
// based on local deformation regime: convergent, divergent, or neutral.
// This reflects physical differences in weakening due to fluids, melt, or
// the absence of both. The dynamic angle is selected accordingly.

const double convergence_indicator = -trace(strain_rate);

double target_angle;

if (convergence_indicator > convergence_threshold)
target_angle = dynamic_angles_of_internal_friction_for_convergence[volume_fraction_index];
else if (convergence_indicator < -divergence_threshold)
target_angle = dynamic_angles_of_internal_friction_for_divergence[volume_fraction_index];
else
target_angle = dynamic_angles_of_internal_friction[volume_fraction_index];

const double mu = std::tan(target_angle)
+ (std::tan(static_friction_angle) - std::tan(target_angle))
/ (1. + std::pow((current_edot_ii / dynamic_characteristic_strain_rate),
dynamic_friction_smoothness_exponent));
const double dynamic_friction_angle = std::atan(mu);
Assert((mu < 1) && (0 < dynamic_friction_angle) && (dynamic_friction_angle <= 1.6), ExcMessage(
"The friction coefficient should be larger than zero and smaller than 1. "
"The friction angle should be smaller than 1.6 rad."));
return dynamic_friction_angle;
}
case function:
{
// Use a given function input per composition to get the friction angle
Expand Down Expand Up @@ -113,24 +154,22 @@ namespace aspect
return static_friction_angle;
}



template <int dim>
FrictionMechanism
FrictionModels<dim>::
get_friction_mechanism() const
{
return friction_mechanism;
}



template <int dim>
void
FrictionModels<dim>::declare_parameters (ParameterHandler &prm)
{
prm.declare_entry ("Friction mechanism", "none",
Patterns::Selection("none|dynamic friction|function"),
Patterns::Selection("none|dynamic friction|differential dynamic friction|function"),
"Whether to make the friction angle dependent on strain rate or not. This rheology "
"is intended to be used together with the visco-plastic rheology model."
"\n\n"
Expand Down Expand Up @@ -172,17 +211,33 @@ namespace aspect
"those corresponding to chemical compositions. "
"Dynamic angles of friction are used as the current friction angle when the effective "
"strain rate is well above the 'dynamic characteristic strain rate'. "
"Units: \\si{\\degree}.");
"Units: \\si{\\degree}.");

prm.declare_entry ("Dynamic friction smoothness exponent", "1",
Patterns::Double (0),
Patterns::List(Patterns::Double(0)),
"An exponential factor in the equation for the calculation of the friction angle when a "
"static and a dynamic angle of internal friction are specified. A factor of 1 returns the equation "
"to Equation (13) in \\cite{van_dinther_seismic_2013}. A factor between 0 and 1 makes the "
"curve of the friction angle vs. the strain rate smoother, while a factor $>$ 1 makes "
"the change between static and dynamic friction angle more steplike. "
"Units: none.");

prm.declare_entry ("Convergence threshold", "1e-15",
Patterns::Double(),
"Minimum convergence rate (negative strain rate trace) to trigger convergent friction angle.");

prm.declare_entry ("Divergence threshold", "1e-15",
Patterns::Double(),
"Minimum divergence rate to trigger divergent friction angle.");

prm.declare_entry ("Dynamic angles of internal friction for convergent flow", "2",
Patterns::List(Patterns::Double(0)),
"Optional dynamic friction angles for convergent flow.");

prm.declare_entry ("Dynamic angles of internal friction for divergent flow", "4",
Patterns::List(Patterns::Double(0)),
"Optional dynamic friction angles for divergent flow.");

/**
* If friction is specified as a function input.
*/
Expand Down Expand Up @@ -210,8 +265,6 @@ namespace aspect
prm.leave_subsection();
}



template <int dim>
void
FrictionModels<dim>::parse_parameters (ParameterHandler &prm)
Expand All @@ -221,6 +274,8 @@ namespace aspect
friction_mechanism = static_friction;
else if (prm.get ("Friction mechanism") == "dynamic friction")
friction_mechanism = dynamic_friction;
else if (prm.get ("Friction mechanism") == "differential dynamic friction")
friction_mechanism = differential_dynamic_friction;
else if (prm.get ("Friction mechanism") == "function")
friction_mechanism = function;
else
Expand All @@ -229,6 +284,11 @@ namespace aspect
// Dynamic friction parameters
dynamic_characteristic_strain_rate = prm.get_double("Dynamic characteristic strain rate");

// Differential dynamic friction parameters
convergence_threshold = prm.get_double("Convergence threshold");
divergence_threshold = prm.get_double("Divergence threshold");


// Retrieve the list of composition names
std::vector<std::string> compositional_field_names = this->introspection().get_composition_names();

Expand All @@ -240,23 +300,47 @@ namespace aspect
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, "Dynamic angles of internal friction");
options.list_of_allowed_keys = compositional_field_names;
Utilities::MapParsing::Options options_default(chemical_field_names, "Dynamic angles of internal friction");
options_default.list_of_allowed_keys = compositional_field_names;

dynamic_angles_of_internal_friction = Utilities::MapParsing::parse_map_to_double_array (prm.get("Dynamic angles of internal friction"),
options);
options_default);

// --- Convergent dynamic angles ---
Utilities::MapParsing::Options options_convergence(chemical_field_names, "Dynamic angles of internal friction for convergent flow");
options_convergence.list_of_allowed_keys = compositional_field_names;

dynamic_angles_of_internal_friction_for_convergence = Utilities::MapParsing::parse_map_to_double_array (prm.get("Dynamic angles of internal friction for convergent flow"),
options_convergence);

// --- Divergent dynamic angles ---
Utilities::MapParsing::Options options_divergence(chemical_field_names, "Dynamic angles of internal friction for divergent flow");
options_divergence.list_of_allowed_keys = compositional_field_names;

dynamic_angles_of_internal_friction_for_divergence = Utilities::MapParsing::parse_map_to_double_array (prm.get("Dynamic angles of internal friction for divergent flow"),
options_divergence);


// Convert angles from degrees to radians
for (double &angle : dynamic_angles_of_internal_friction)
{
AssertThrow(angle <= 90,
ExcMessage("Dynamic angles of friction must be <= 90 degrees"));
angle *= constants::degree_to_radians;
}
auto convert_angles_to_radians = [](std::vector<double> &angles, const std::string &description)
{
for (double &angle : angles)
{
AssertThrow(angle <= 90,
ExcMessage(description + " must be <= 90 degrees"));
angle *= constants::degree_to_radians;
}
};

convert_angles_to_radians(dynamic_angles_of_internal_friction, "Dynamic angles of internal friction");
convert_angles_to_radians(dynamic_angles_of_internal_friction_for_convergence, "Dynamic angles of internal friction for convergent flow");
convert_angles_to_radians(dynamic_angles_of_internal_friction_for_divergence, "Dynamic angles of internal friction for divergent flow");


dynamic_friction_smoothness_exponent = prm.get_double("Dynamic friction smoothness exponent");



// Get the number of fields for composition-dependent material properties
// including the background field.
// TODO Make sure functions only have to be specified per chemical composition,
Expand Down
3 changes: 2 additions & 1 deletion source/material_model/rheology/visco_plastic.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -335,7 +335,8 @@ namespace aspect
output_parameters.drucker_prager_parameters[j].angle_internal_friction = friction_models.compute_friction_angle(effective_edot_ii,
j,
output_parameters.drucker_prager_parameters[j].angle_internal_friction,
in.position[i]);
in.position[i],
in.strain_rate[i]);

// Step 5: plastic yielding

Expand Down