Changing to use viscous viscosity for tidal heating rate calculation when using Visco Plastic material model by hyunseong96 · Pull Request #7168 · geodynamics/aspect

hyunseong96 · 2026-07-28T09:16:24Z

Hiya, I create a PR to change viscosity used to calculate tidal heating rate in visco_plastic.cc.
The heating rate equation follows eq 12 of Tobie et al. (2003) (https://doi.org/10.1029/2003JE002099).

where $\eta$ is viscosity, $\bar{\dot{\epsilon}^2}$ is the time-averaged square of tidal strain rate, $\omega$ is tidal frequency, and $\mu_E$ is the elastic shear modulus.

Currently, ASPECT uses output viscosity. It means, when elasticity is enabled, effective elasticity is used. However, tidal heating rate equation is already a simplified viscoelastic shear heating equation and it is aimed for the period of orbiting. It means the process of making VE viscosity is used twice; one time for making VE viscosity of convective time scale and the other time for making VE viscosity of orbiting time scale. Also, it is likely that elastic stress that is accumulated over convective time scale does not affect day-scale deformation of tidal heating. Therefore, using effective viscosity for the tidal heating rate calculation only implies lowered viscosity causing undesirable larger tidal heating.

@Geoniette and I had a discussion and decided to change the viscosity used for tidal heating as viscous viscosity.
In visco_plastic.cc of rheology, I made composition_viscous_viscosities that saves non_yielding_viscosity right before viscoelasticity is included to non_yielding_viscosity and that is truncated by maximum and minimum viscosity settings. Then, composition_viscous_viscosities is saved to additional output which will be used for tidal heating calculation in visco_plastic.cc of material_model. I put if-statement to use this viscous_viscosities when elasticity is enabled in tidal_heating.cc.

The changes follow as the figure below. The result is a simple setting of ice with linear temperature increase by depth at time=0. First row shows the tidal heating rate, and second row shows the viscosity. Two left columns show the viscosities and heating rates based on current implementation under viscoelasticity and viscosity. Two right columns shows the heating rate and viscosity based on new implementation under viscoelasticity with two different elastic shear modulus. Comparing first and second columns, you can see that tidal heating rate is larger when elasticity is included as mentioned. Third column is the new implementation of same viscoelasticity setting with the first column and shows same tidal heating rate with second column of viscous deformation setting. Third and fourth columns show tidal heating rate under different elastic shear modulus (1e9 and 1e7 Pa). These show that tidal heating rate is dependent on viscous viscosity, regardless of elasticity setting!

Cheers,
Hyunseong

bangerth · 2026-07-28T14:33:43Z

You have some merge conflicts. Can you resolve those?

hyunseong96 · 2026-07-28T15:38:34Z

@bangerth Yes I have check the conflicts and resolved them! There was a recent merge on visco_plastic.cc while I was working on few days older version than that.

gassmoeller

I think this PR needs some structural changes. I am ok with the feature you are aiming for, but the current PR does not follow our usual separation of concerns, and I dont think it is working as you intended. Take a look at my comments and implement them as far as you get. Then find one of us to talk about it.

gassmoeller · 2026-07-28T15:36:18Z

+        /**
+         * Returns true if the tidal heating plugin is found in the
+         * list of active heating models.
+         */
+        bool
+        tidal_heating_enabled() const;


please dont introduce this function, the other functions of this type should be deprecated and not used any longer. I will point out below how to do this differently.

gassmoeller · 2026-07-28T15:37:48Z

+         * outputs.
+         */
+        std::unique_ptr<HeatingModel::TidalHeating<dim>> tidal_heating;
+


This is most likely not what you want. You want to get access to the existing heating model plugin, not create your own copy of a heating model inside the material model. I will point out below how you get access without creating a copy.

gassmoeller · 2026-07-28T15:39:14Z

+    template <int dim>
+    bool
+    Manager<dim>::tidal_heating_enabled() const
+    {
+      return this->template has_matching_active_plugin<HeatingModel::TidalHeating<dim>>() ;
+    }


Oh, so you know already how to do this. No need to create this function then. Just call this->get_heating_model_manager().template has_matching_active_plugin<HeatingModel::TidalHeating<dim>>() in the place where you call this function at the moment.

gassmoeller · 2026-07-28T15:41:51Z

+          double viscous_viscosities = material_model_outputs.viscosities[q];
+
+          if (this->get_parameters().enable_elasticity)
+            {
+              const std::shared_ptr<const HeatingModel::ViscousAdditionalOutputs<dim>> viscous_out =
+                material_model_outputs.template get_additional_output_object<HeatingModel::ViscousAdditionalOutputs<dim>>();
+              if (viscous_out != nullptr)
+                {
+                  viscous_viscosities = viscous_out->viscous_viscosity[q];
+                }
+            }
+
+
+          heating_model_outputs.heating_source_terms[q] = 2. * viscous_viscosities * local_tidal_strain_rate * local_tidal_strain_rate
+                                                          / ( 1. + ( (tidal_frequency * viscous_viscosities) / elastic_shear_modulus ) * ( (tidal_frequency * viscous_viscosities) / elastic_shear_modulus ) );


viscous_viscosity always sounds silly to me. And I think instead you could write

// determine the viscosity, if elasticity is enabled, only consider viscous deformation double viscosity = material_model_outputs.viscosities[q]; if (...) ... viscosity = ....

gassmoeller · 2026-07-28T15:44:16Z

    }


-


did you delete this line on purpose? please reintroduce it

gassmoeller · 2026-07-28T15:45:09Z


-
-            // Step 1e: multiply the viscosity by a constant (default value is 1)
+            // Step 1f: multiply the viscosity by a constant (default value is 1)


this is a correct fix, but it is unrelated to this PR. Please open a separate PR for it (it will be quick to merge).

gassmoeller · 2026-07-28T15:57:40Z

+              // Fill the additional viscous outputs that will be used for tidal heating.
+              if (in.requests_property(MaterialProperties::additional_outputs) && this->get_parameters().enable_elasticity)
+                tidal_heating->fill_additional_viscous_outputs(i, volume_fractions, out, isostrain_viscosities, rheology->viscosity_averaging);
+


but at the moment you do not know if tidal heating is actually used in this model. This compiles at the moment because you always create a copy of tidal heating inside visco plastic, but this is not the correct approach. What you really want to do conceptually is:

if (in.requests_property(MaterialProperties::additional_outputs) && viscous_outputs exist) .. fill additional viscous outputs

The conditional statement would look like this:

if (in.requests_property(MaterialProperties::additional_outputs) && out.template has_additional_output_object<MaterialModel::ViscousAdditionalOutputs<dim>>()) { fill_additional_viscous_outputs(i, volume_fractions, out, isostrain_viscosities, rheology->viscosity_averaging); } And the fill additional viscous outputs function should live inside this class (it is a feature that this class offers, it is not a feature of tidal heating to know how to compute viscous viscosity). There is nothing that makes viscous outputs specific to the tidal heating model, it just happens to be the single plugin that needs this output at the moment. But that may change in the future (and then we do not want to depend on tidal heating if it is not even active in that model). So move everything possible into this class.

gassmoeller · 2026-07-28T15:59:45Z

+    TidalHeating<dim>::create_additional_material_model_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
+    {
+      if (out.template has_additional_output_object<ViscousAdditionalOutputs<dim>>() == false)
+        {
+          const unsigned int n_points = out.n_evaluation_points();
+          out.additional_outputs.push_back(
+            std::make_unique<ViscousAdditionalOutputs<dim>> (n_points));
+        }
+    }


This is the only function of the ones below that should stay in this file and class.

gassmoeller · 2026-07-28T16:00:05Z

+    template <int dim>
+    ViscousAdditionalOutputs<dim>::ViscousAdditionalOutputs(const unsigned int n_points)
+      : MaterialModel::NamedAdditionalMaterialOutputs<dim>({"viscous_viscosity"}),
+    viscous_viscosity(n_points, std::numeric_limits<double>::max())
+    {}
+
+    template <int dim>
+    std::vector<double>
+    ViscousAdditionalOutputs<dim>::get_nth_output(const unsigned int idx) const
+    {
+      (void) idx;
+      AssertIndexRange (idx, 1);
+
+      return viscous_viscosity;
+    }
+
+    template <int dim>
+    void
+    TidalHeating<dim>::fill_additional_viscous_outputs (const unsigned int i,
+                                                        const std::vector<double> &volume_fractions,
+                                                        MaterialModel::MaterialModelOutputs<dim> &out,
+                                                        const MaterialModel::IsostrainViscosities &isostrain_viscosities,
+                                                        const MaterialModel::MaterialUtilities::CompositionalAveragingOperation &average_type) const
+    {
+      const std::shared_ptr<ViscousAdditionalOutputs<dim>> viscous_out =
+        out.template get_additional_output_object<ViscousAdditionalOutputs<dim>>();
+      if (viscous_out != nullptr)
+        {
+          viscous_out->viscous_viscosity[i] =
+            MaterialModel::MaterialUtilities::average_value(volume_fractions,
+                                                            isostrain_viscosities.composition_viscous_viscosities,
+                                                            average_type);
+        }
+    }


All of this should move over to visco_plastic.

gassmoeller · 2026-07-28T16:02:54Z

        output_parameters.composition_yielding.resize(volume_fractions.size(), false);
        output_parameters.composition_viscosities.resize(volume_fractions.size(), numbers::signaling_nan<double>());
+        output_parameters.composition_viscous_viscosities.resize(volume_fractions.size(), numbers::signaling_nan<double>());
        output_parameters.drucker_prager_parameters.resize(volume_fractions.size());
        output_parameters.dilation_lhs_terms.resize(volume_fractions.size(), numbers::signaling_nan<double>());
        output_parameters.dilation_rhs_terms.resize(volume_fractions.size(), numbers::signaling_nan<double>());


we need to find a way to avoid accumulating more and more vectors in this structure. I commented on that in #7133. Whatever the solution to that PR will also need to be applied here.

hyunseong96 · 2026-07-29T08:54:34Z

Hi @gassmoeller ,
I have moved functions related to additional output except create function and removed unused headers and lines that I put.

hyunseong96 · 2026-07-29T19:11:59Z

Hi @gassmoeller ,
I have moved functions related to generating additional outputs to material_model/additional_outputs folder.
Also I changed the names of parameters and function names as #7133 has similar function names.
Please let me know for further changes :)
Thank you very much!

gassmoeller

Yes, the direction of this PR is correct. I have some comments on file names, and I still have not thought of a good solution for the IsostrainViscosities structure, but I need to think about that when I am more awake. Please for now implement my suggestions, they should get you to a working state, which is mostly complete.

gassmoeller · 2026-07-29T19:35:24Z

keep the name of the file consistent with the name of the class, not the purpose you are using it for. This class can be used by other users for other purposes.

I would suggest renaming to viscosity_without_elasticity

gassmoeller · 2026-07-29T19:35:51Z

+#ifndef _aspect_material_model_additional_outputs_viscosity_for_tidal_heating_h
+#define _aspect_material_model_additional_outputs_viscosity_for_tidal_heating_h


viscosity_without_elasticity in both lines

gassmoeller · 2026-07-29T19:36:05Z

+#ifndef _aspect_material_model_additional_outputs_viscosity_for_tidal_heating_h
+#define _aspect_material_model_additional_outputs_viscosity_for_tidal_heating_h
+
+#include <aspect/simulator_access.h>


do you need simulator_access?

gassmoeller · 2026-07-29T19:36:21Z

+
+#include <aspect/simulator_access.h>
+#include <aspect/material_model/interface.h>
+#include <aspect/material_model/visco_plastic.h>


you definitely do not need visco_plastic here

gassmoeller · 2026-07-29T19:38:21Z

 #include <aspect/heating_model/tidal_heating.h>

+#include <aspect/material_model/additional_outputs/viscosity_for_tidal_heating.h>
+#include <aspect/material_model/visco_plastic.h>


i dont think you need the visco_plastic header here

gassmoeller · 2026-07-29T19:41:08Z

+       * viscosity = material_model_outputs.viscosities when elasticity is not enabled. If enabled, material_model_outputs.viscosities_without_elasticity.
+       * Because tidal heating equation is already a form of viscoelastic shear heating while orbiting,


Please rewrite this so that it explains viscosity. I.e. something like
viscosity should only contain the viscous part of the deformation rate, i.e. not plastic or elastic contributions to the effective viscosity. This is because tidal heating is a form of viscous dissipation that only occurs due to viscous deformation. (or something similar, maybe I misunderstood the justification)

gassmoeller · 2026-07-29T19:41:49Z

same as above, rename the file to align with the class name

gassmoeller · 2026-07-29T19:44:03Z

+
+    template <int dim>
+    void
+    ViscoPlastic<dim>::fill_viscosity_for_tidal_heating_outputs (const unsigned int i,


this function should be called fill_viscosity_without_elasticity_outputs.

gassmoeller · 2026-07-29T19:45:32Z

+Changed: Now tidal heating module only uses viscous viscosity regardless of enabling elasticity in material model Visco Plastic.
+Because the tidal heating rate equation is already a form of viscoelastic shear heating while orbiting,
+inserting viscoelastic viscosity again to the shear heating equation is incorrect.


Suggested change

Changed: Now tidal heating module only uses viscous viscosity regardless of enabling elasticity in material model Visco Plastic.

Because the tidal heating rate equation is already a form of viscoelastic shear heating while orbiting,

inserting viscoelastic viscosity again to the shear heating equation is incorrect.

Changed: Now the tidal heating module only uses the viscosity without elastic contributions from the material model visco plastic.

hyunseong96 force-pushed the tidal_heating_viscosity branch from 882ba7c to 24771d3 Compare July 28, 2026 14:36

gassmoeller reviewed Jul 28, 2026

View reviewed changes

hyunseong96 force-pushed the tidal_heating_viscosity branch 3 times, most recently from da51363 to 6dae6e2 Compare July 29, 2026 07:45

hyunseong96 force-pushed the tidal_heating_viscosity branch from 6dae6e2 to dcf44dd Compare July 29, 2026 09:57

gassmoeller added ready to merge ready to test and removed ready to merge labels Jul 29, 2026

hyunseong96 force-pushed the tidal_heating_viscosity branch from dcf44dd to 5f36abf Compare July 29, 2026 19:00

Getting viscosity of viscous rheology for tidal heating calculation

43d3330

hyunseong96 force-pushed the tidal_heating_viscosity branch from 5f36abf to fe86f15 Compare July 29, 2026 19:11

hyunseong96 force-pushed the tidal_heating_viscosity branch from fe86f15 to 4390fd1 Compare July 29, 2026 19:29

gassmoeller reviewed Jul 29, 2026

View reviewed changes

Moving functions to material_model/additional_outputs

87e3814

hyunseong96 force-pushed the tidal_heating_viscosity branch from e5b0834 to 87e3814 Compare July 30, 2026 05:28



		// Step 1e: multiply the viscosity by a constant (default value is 1)
		// Step 1f: multiply the viscosity by a constant (default value is 1)

		#ifndef _aspect_material_model_additional_outputs_viscosity_for_tidal_heating_h
		#define _aspect_material_model_additional_outputs_viscosity_for_tidal_heating_h

		* viscosity = material_model_outputs.viscosities when elasticity is not enabled. If enabled, material_model_outputs.viscosities_without_elasticity.
		* Because tidal heating equation is already a form of viscoelastic shear heating while orbiting,

Uh oh!

Conversation

hyunseong96 commented Jul 28, 2026

Uh oh!

bangerth commented Jul 28, 2026

Uh oh!

hyunseong96 commented Jul 28, 2026

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hyunseong96 commented Jul 29, 2026

Uh oh!

hyunseong96 commented Jul 29, 2026

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants