diff --git a/doc/modules/changes/20260728_Minerallo b/doc/modules/changes/20260728_Minerallo
new file mode 100644
index 00000000000..7d15f7a5794
--- /dev/null
+++ b/doc/modules/changes/20260728_Minerallo
@@ -0,0 +1,6 @@
+Fixed: When mesh deformation is enabled, the adiabatic pressure used by the
+visco-plastic material model for plastic yielding and viscous creep is now
+evaluated using depth below the current deformed surface. This prevents surface
+depressions such as trenches from becoming artificially strong.
+
+(Michaël Pons, 2026/07/28)
diff --git a/doc/sphinx/parameters/Material_20model.md b/doc/sphinx/parameters/Material_20model.md
index 5bdb25f020a..127b6b3be83 100644
--- a/doc/sphinx/parameters/Material_20model.md
+++ b/doc/sphinx/parameters/Material_20model.md
@@ -5547,7 +5547,7 @@ If a compositional field named ’noninitial\_plastic\_strain’ is incl
**Pattern:** [Bool]
-**Documentation:** Whether to use the adiabatic pressure instead of the full pressure (default) when calculating viscous creep. This may be helpful in models where the full pressure has an unusually large negative value arising from large negative dynamic pressure, resulting in solver convergence issue and in some cases a viscosity of zero.
+**Documentation:** Whether to use the adiabatic pressure instead of the full pressure (default) when calculating viscous creep. This may be helpful in models where the full pressure has an unusually large negative value arising from large negative dynamic pressure, resulting in solver convergence issue and in some cases a viscosity of zero. With mesh deformation, the adiabatic profile is evaluated using depth below the current deformed surface.
::::
::::{dropdown} __Parameter:__ {ref}`Use adiabatic pressure in plasticity`
@@ -5556,7 +5556,7 @@ If a compositional field named ’noninitial\_plastic\_strain’ is incl
**Pattern:** [Bool]
-**Documentation:** Whether to use the adiabatic pressure instead of the full pressure when calculating plastic yield stress. This may be helpful in models where the full pressure has unusually large variations, resulting in solver convergence issues. Be aware that this setting will change the plastic shear band angle.
+**Documentation:** Whether to use the adiabatic pressure instead of the full pressure when calculating plastic yield stress. This may be helpful in models where the full pressure has unusually large variations, resulting in solver convergence issues. With mesh deformation, the adiabatic profile is evaluated using depth below the current deformed surface, so surface depressions such as oceanic trenches remain at the adiabatic surface pressure. Be aware that this setting will also change the plastic shear band angle.
::::
::::{dropdown} __Parameter:__ {ref}`Use dominant phase for viscosity`
diff --git a/include/aspect/adiabatic_conditions/interface.h b/include/aspect/adiabatic_conditions/interface.h
index 3d494864bbe..bdedf8bf8f9 100644
--- a/include/aspect/adiabatic_conditions/interface.h
+++ b/include/aspect/adiabatic_conditions/interface.h
@@ -76,6 +76,14 @@ namespace aspect
virtual
double pressure (const Point &p) const = 0;
+ /**
+ * Return the adiabatic pressure at the representative point
+ * corresponding to @p depth. The depth must lie between zero and the
+ * maximal depth of the geometry model.
+ */
+ double
+ pressure_at_depth (const double depth) const;
+
/**
* Return the reference_density at a given point of the domain.
*/
diff --git a/include/aspect/material_model/rheology/visco_plastic.h b/include/aspect/material_model/rheology/visco_plastic.h
index 40b8a4cf6a7..3530f74218e 100644
--- a/include/aspect/material_model/rheology/visco_plastic.h
+++ b/include/aspect/material_model/rheology/visco_plastic.h
@@ -204,11 +204,23 @@ namespace aspect
IsostrainViscosities
calculate_isostrain_viscosities ( const MaterialModel::MaterialModelInputs &in,
const unsigned int i,
+ const double current_surface_adiabatic_pressure,
const std::vector &volume_fractions,
const std::vector &phase_function_values = std::vector(),
const std::vector &n_phase_transitions_per_composition =
std::vector()) const;
+ /**
+ * Return adiabatic pressures evaluated using depth below the current
+ * surface at all evaluation points. These values are used by the
+ * pressure-dependent parts of the rheology for which adiabatic
+ * pressure is enabled. A NaN entry indicates that the existing
+ * pressure selection in the rheology should be used.
+ */
+ std::vector
+ compute_current_surface_adiabatic_pressures(
+ const MaterialModel::MaterialModelInputs &in) const;
+
/**
* A function that fills the viscosity derivatives in the
* MaterialModelOutputs object that is handed over, if they exist.
@@ -221,6 +233,7 @@ namespace aspect
void compute_viscosity_derivatives(const unsigned int point_index,
const std::vector &volume_fractions,
const IsostrainViscosities &isostrain_values,
+ const double current_surface_adiabatic_pressure,
const MaterialModel::MaterialModelInputs &in,
MaterialModel::MaterialModelOutputs &out,
const std::vector &phase_function_values = std::vector(),
@@ -268,6 +281,7 @@ namespace aspect
void fill_plastic_outputs(const unsigned int point_index,
const std::vector &volume_fractions,
const bool plastic_yielding,
+ const double current_surface_adiabatic_pressure,
const MaterialModel::MaterialModelInputs &in,
MaterialModel::MaterialModelOutputs &out,
const IsostrainViscosities &isostrain_viscosities) const;
diff --git a/source/adiabatic_conditions/interface.cc b/source/adiabatic_conditions/interface.cc
index 8097f2f2e91..e657bd66ae9 100644
--- a/source/adiabatic_conditions/interface.cc
+++ b/source/adiabatic_conditions/interface.cc
@@ -30,6 +30,15 @@ namespace aspect
{
namespace AdiabaticConditions
{
+ template
+ double
+ Interface::pressure_at_depth(const double depth) const
+ {
+ return pressure(this->get_geometry_model().representative_point(depth));
+ }
+
+
+
template
void Interface::get_adiabatic_temperature_profile(std::vector &values) const
{
diff --git a/source/material_model/rheology/visco_plastic.cc b/source/material_model/rheology/visco_plastic.cc
index 2e88662367e..b8e8c10deca 100644
--- a/source/material_model/rheology/visco_plastic.cc
+++ b/source/material_model/rheology/visco_plastic.cc
@@ -19,6 +19,7 @@
*/
#include
+#include
#include
#include
#include
@@ -158,6 +159,7 @@ namespace aspect
ViscoPlastic::
calculate_isostrain_viscosities (const MaterialModel::MaterialModelInputs &in,
const unsigned int i,
+ const double current_surface_adiabatic_pressure,
const std::vector &volume_fractions,
const std::vector &phase_function_values,
const std::vector &n_phase_transitions_per_composition) const
@@ -250,7 +252,12 @@ namespace aspect
double pressure_for_creep = in.pressure[i];
if (use_adiabatic_pressure_in_creep)
- pressure_for_creep = this->get_adiabatic_conditions().pressure(in.position[i]);
+ pressure_for_creep =
+ (std::isnan(current_surface_adiabatic_pressure)
+ ?
+ this->get_adiabatic_conditions().pressure(in.position[i])
+ :
+ current_surface_adiabatic_pressure);
const double viscosity_diffusion
= (viscous_flow_law != dislocation
@@ -429,9 +436,13 @@ namespace aspect
// than the lithostatic pressure.
double pressure_for_plasticity = in.pressure[i];
-
if (use_adiabatic_pressure_in_plasticity)
- pressure_for_plasticity = this->get_adiabatic_conditions().pressure(in.position[i]);
+ pressure_for_plasticity =
+ (std::isnan(current_surface_adiabatic_pressure)
+ ?
+ this->get_adiabatic_conditions().pressure(in.position[i])
+ :
+ current_surface_adiabatic_pressure);
if (allow_negative_pressures_in_plasticity == false)
pressure_for_plasticity = std::max(pressure_for_plasticity,0.0);
@@ -514,12 +525,48 @@ namespace aspect
+ template
+ std::vector
+ ViscoPlastic::compute_current_surface_adiabatic_pressures(
+ const MaterialModel::MaterialModelInputs &in) const
+ {
+ std::vector pressures(
+ in.n_evaluation_points(),
+ std::numeric_limits::quiet_NaN());
+
+ if ((use_adiabatic_pressure_in_plasticity == false
+ && use_adiabatic_pressure_in_creep == false)
+ || this->get_parameters().mesh_deformation_enabled == false
+ || in.current_cell.state() != IteratorState::valid)
+ return pressures;
+
+ const std::vector depths =
+ this->get_mesh_deformation_handler().depth_below_current_surface(
+ in.current_cell, in.position);
+ const double maximal_depth = this->get_geometry_model().maximal_depth();
+
+ for (unsigned int q = 0; q < in.n_evaluation_points(); ++q)
+ {
+ // A deformed domain can extend beyond the reference geometry.
+ // representative_point() is only defined inside its depth range.
+ const double bounded_depth =
+ std::min(std::max(depths[q], 0.0), maximal_depth);
+ pressures[q] =
+ this->get_adiabatic_conditions().pressure_at_depth(bounded_depth);
+ }
+
+ return pressures;
+ }
+
+
+
template
void
ViscoPlastic::
compute_viscosity_derivatives(const unsigned int i,
const std::vector &volume_fractions,
const IsostrainViscosities ¤t_isostrain_values,
+ const double current_surface_adiabatic_pressure,
const MaterialModel::MaterialModelInputs &in,
MaterialModel::MaterialModelOutputs &out,
const std::vector &phase_function_values,
@@ -565,7 +612,9 @@ namespace aspect
in_derivatives.strain_rate[i] = forward_strain_rate;
const IsostrainViscosities forward_isostrain_values =
- calculate_isostrain_viscosities(in_derivatives, i, volume_fractions,
+ calculate_isostrain_viscosities(in_derivatives, i,
+ current_surface_adiabatic_pressure,
+ volume_fractions,
phase_function_values, n_phase_transitions_per_composition);
// For each composition of the independent component, compute the derivative.
@@ -603,7 +652,9 @@ namespace aspect
in_derivatives.strain_rate[i] = current_strain_rate;
const IsostrainViscosities forward_isostrain_values =
- calculate_isostrain_viscosities(in_derivatives, i, volume_fractions,
+ calculate_isostrain_viscosities(in_derivatives, i,
+ current_surface_adiabatic_pressure,
+ volume_fractions,
phase_function_values, n_phase_transitions_per_composition);
for (unsigned int composition_index = 0; composition_index < n_compositions; ++composition_index)
@@ -769,15 +820,21 @@ namespace aspect
"This may be helpful in models where the "
"full pressure has an unusually large negative value arising from "
"large negative dynamic pressure, resulting in solver convergence "
- "issue and in some cases a viscosity of zero.");
+ "issue and in some cases a viscosity of zero. With mesh deformation, "
+ "the adiabatic profile is evaluated using depth below the current "
+ "deformed surface.");
prm.declare_entry ("Use adiabatic pressure in plasticity", "false",
Patterns::Bool (),
"Whether to use the adiabatic pressure instead of the full "
"pressure when calculating plastic yield stress. "
"This may be helpful in models where the "
"full pressure has unusually large variations, resulting "
- "in solver convergence issues. Be aware that this setting "
- "will change the plastic shear band angle.");
+ "in solver convergence issues. With mesh deformation, the "
+ "adiabatic profile is evaluated using depth below the current "
+ "deformed surface, so surface depressions such as oceanic "
+ "trenches remain at the adiabatic surface pressure. Be aware "
+ "that this setting will also change the "
+ "plastic shear band angle.");
// Diffusion creep parameters
Rheology::DiffusionCreep::declare_parameters(prm);
@@ -1014,6 +1071,7 @@ namespace aspect
fill_plastic_outputs(const unsigned int i,
const std::vector &volume_fractions,
const bool plastic_yielding,
+ const double current_surface_adiabatic_pressure,
const MaterialModel::MaterialModelInputs &in,
MaterialModel::MaterialModelOutputs &out,
const IsostrainViscosities &isostrain_viscosities) const
@@ -1032,9 +1090,13 @@ namespace aspect
plastic_out->yielding[i] = plastic_yielding ? 1 : 0;
double pressure_for_plasticity = in.pressure[i];
-
if (use_adiabatic_pressure_in_plasticity)
- pressure_for_plasticity = this->get_adiabatic_conditions().pressure(in.position[i]);
+ pressure_for_plasticity =
+ (std::isnan(current_surface_adiabatic_pressure)
+ ?
+ this->get_adiabatic_conditions().pressure(in.position[i])
+ :
+ current_surface_adiabatic_pressure);
if (allow_negative_pressures_in_plasticity == false)
pressure_for_plasticity = std::max(pressure_for_plasticity, 0.0);
diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc
index 7959eee6224..fcfb50825ff 100644
--- a/source/material_model/visco_plastic.cc
+++ b/source/material_model/visco_plastic.cc
@@ -88,7 +88,14 @@ namespace aspect
/* The following returns whether or not the material is plastically yielding
* as documented in evaluate.
*/
- const IsostrainViscosities isostrain_viscosities = rheology->calculate_isostrain_viscosities(in, 0, volume_fractions, phase_function_values, phase_function.n_phase_transitions_for_each_composition());
+ const std::vector current_surface_adiabatic_pressures =
+ rheology->compute_current_surface_adiabatic_pressures(in);
+ const IsostrainViscosities isostrain_viscosities =
+ rheology->calculate_isostrain_viscosities(in, 0,
+ current_surface_adiabatic_pressures[0],
+ volume_fractions,
+ phase_function_values,
+ phase_function.n_phase_transitions_for_each_composition());
std::vector::const_iterator max_composition = std::max_element(volume_fractions.begin(), volume_fractions.end());
const bool plastic_yielding = isostrain_viscosities.composition_yielding[std::distance(volume_fractions.begin(), max_composition)];
@@ -116,6 +123,8 @@ namespace aspect
std::vector phase_function_discrete_values = (use_dominant_phase_for_viscosity?
std::vector(phase_function_discrete->n_phase_transitions(), 0.0): std::vector());
+ const std::vector current_surface_adiabatic_pressures =
+ rheology->compute_current_surface_adiabatic_pressures(in);
// Loop through all requested points
for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
@@ -214,12 +223,20 @@ namespace aspect
phase_function_discrete_values[j] = phase_function_discrete->compute_value(phase_inputs);
}
isostrain_viscosities =
- rheology->calculate_isostrain_viscosities(in, i, volume_fractions, phase_function_discrete_values, phase_function_discrete->n_phase_transitions_for_each_chemical_composition());
+ rheology->calculate_isostrain_viscosities(in, i,
+ current_surface_adiabatic_pressures[i],
+ volume_fractions,
+ phase_function_discrete_values,
+ phase_function_discrete->n_phase_transitions_for_each_chemical_composition());
}
else
{
isostrain_viscosities =
- rheology->calculate_isostrain_viscosities(in, i, volume_fractions, phase_function_values, n_phase_transitions_for_each_chemical_composition);
+ rheology->calculate_isostrain_viscosities(in, i,
+ current_surface_adiabatic_pressures[i],
+ volume_fractions,
+ phase_function_values,
+ n_phase_transitions_for_each_chemical_composition);
}
// The isostrain condition implies that the viscosity averaging should be arithmetic (see above).
@@ -241,6 +258,7 @@ namespace aspect
rheology->compute_viscosity_derivatives(i, volume_fractions,
isostrain_viscosities,
+ current_surface_adiabatic_pressures[i],
in, out, phase_function_values,
n_phase_transitions_for_each_chemical_composition);
}
@@ -277,7 +295,9 @@ namespace aspect
// has been called.
if (in.requests_property(MaterialProperties::additional_outputs))
{
- rheology->fill_plastic_outputs(i, volume_fractions, plastic_yielding, in, out, isostrain_viscosities);
+ rheology->fill_plastic_outputs(i, volume_fractions, plastic_yielding,
+ current_surface_adiabatic_pressures[i],
+ in, out, isostrain_viscosities);
rheology->fill_viscosity_outputs(i, volume_fractions, out, isostrain_viscosities);
}
diff --git a/tests/visco_plastic_adiabatic_pressure_mesh_deformation.prm b/tests/visco_plastic_adiabatic_pressure_mesh_deformation.prm
new file mode 100644
index 00000000000..424afce53e6
--- /dev/null
+++ b/tests/visco_plastic_adiabatic_pressure_mesh_deformation.prm
@@ -0,0 +1,99 @@
+# Verify that adiabatic pressure used by plasticity and viscous creep is
+# evaluated relative to the current surface. The right half of the top boundary
+# moves down by 6 km. Yield stress and diffusion-creep viscosity along both
+# halves of the final surface should remain at their zero-pressure values.
+
+set Dimension = 2
+set End time = 1e4
+set Maximum time step = 1e4
+set Use years instead of seconds = true
+set Pressure normalization = no
+set Surface pressure = 0
+
+subsection Adiabatic conditions model
+ set Model name = function
+
+ subsection Function
+ set Function expression = 1573; 3.3e4*depth; 3300
+ end
+end
+
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X extent = 100e3
+ set Y extent = 50e3
+ set X repetitions = 2
+ set Y repetitions = 1
+ end
+end
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Initial global refinement = 2
+end
+
+subsection Mesh deformation
+ set Mesh deformation boundary indicators = top : free surface & boundary function
+ set Additional tangential mesh velocity boundary indicators = left, right
+
+ subsection Boundary function
+ # Move the right half down by 6 km during the single 10 kyr time step.
+ set Function expression = 0; if(x>50e3, -0.6, 0)
+ end
+end
+
+subsection Boundary temperature model
+ set Fixed temperature boundary indicators = bottom, top
+ set List of model names = initial temperature
+end
+
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = bottom, left, right
+end
+
+subsection Initial temperature model
+ set Model name = function
+
+ subsection Function
+ set Function expression = 1573
+ end
+end
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 0
+ end
+end
+
+subsection Material model
+ set Model name = visco plastic
+
+ subsection Visco Plastic
+ set Reference strain rate = 1e-16
+ set Viscous flow law = composite
+ set Use adiabatic pressure in creep viscosity = true
+ set Use adiabatic pressure in plasticity = true
+ set Angles of internal friction = 10
+ set Cohesions = 1e6
+ end
+end
+
+subsection Formulation
+ set Formulation = Boussinesq approximation
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization
+
+ subsection Visualization
+ set Time between graphical output = 0
+ set Interpolate output = false
+ set List of output variables = named additional outputs, depth, depth including mesh deformation
+ set Output format = gnuplot
+ set Number of grouped files = 1
+ end
+end
diff --git a/tests/visco_plastic_adiabatic_pressure_mesh_deformation.sh b/tests/visco_plastic_adiabatic_pressure_mesh_deformation.sh
new file mode 100755
index 00000000000..b61b1db23c2
--- /dev/null
+++ b/tests/visco_plastic_adiabatic_pressure_mesh_deformation.sh
@@ -0,0 +1,25 @@
+#!/bin/bash
+
+if [ "$1" = "screen-output" ]; then
+ cat > /dev/null
+ echo "Current-surface points (x, y, geometry depth, current depth, pressure used for plasticity in MPa, diffusion viscosity in Pa s):"
+ # Invert the 2D Drucker-Prager yield-stress equation to recover the pressure
+ # used by the material model: yield_stress = C cos(phi) + p sin(phi).
+ awk '
+ BEGIN {
+ pi = atan2(0, -1)
+ cohesion = 1e6
+ friction = 10 * pi / 180
+ }
+ !/^#/ && $14 == 0 {
+ pressure = ($9 - cohesion * cos(friction)) / sin(friction)
+ if (pressure > -500 && pressure < 500)
+ pressure = 0
+ printf "%g %g %.0f %.0f %.3f %.6e\n", $1, $2, $13, $14, pressure / 1e6, $11
+ }' \
+ output-visco_plastic_adiabatic_pressure_mesh_deformation/solution/solution-00001.0000.gnuplot \
+ | sort -n -k1,1 \
+ | uniq
+else
+ cat
+fi
diff --git a/tests/visco_plastic_adiabatic_pressure_mesh_deformation/screen-output b/tests/visco_plastic_adiabatic_pressure_mesh_deformation/screen-output
new file mode 100644
index 00000000000..ba220f1f013
--- /dev/null
+++ b/tests/visco_plastic_adiabatic_pressure_mesh_deformation/screen-output
@@ -0,0 +1,10 @@
+Current-surface points (x, y, geometry depth, current depth, pressure used for plasticity in MPa, diffusion viscosity in Pa s):
+0 50000 0 0 0.000 9.446440e+17
+12500 50000 0 0 0.000 9.446440e+17
+25000 50000 0 0 0.000 9.446440e+17
+37500 50000 0 0 0.000 9.446440e+17
+50000 50000 0 0 0.000 9.446440e+17
+62500 44000 6000 0 0.000 9.446440e+17
+75000 44000 6000 0 0.000 9.446440e+17
+87500 44000 6000 0 0.000 9.446440e+17
+100000 44000 6000 0 0.000 9.446440e+17