m_rheology.py

m_rheology.G(composition)

Returns shear modulus of the material. Can be compositionally-dependent (see the visco-elasticity benchmark in Section 4.2.6).

Parameters:

composition – material composition, e.g. composition = [C_mantle, C_lid, C_plume]

Variables:

G_ice – shear modulus of ice \(G_{ice}\)

Returns:

\[G = G_{ice}\]

or in case the visco-elasticity benchmark

\[G = G_{mantle}^{C_{mantle}} \cdot G_{lid}^{C_{lid}} \cdot G_{plume}^{C_{plume}}\]

m_rheology.cohesion(plastic_strain)

Computes the cohesion depending on the plastic strain (damage) of the material.

Parameters:

• plastic_strain – plastic strain (\(\varepsilon_p\))

• cohesion_strong – lower cut-off value (\(\sigma_Y^{min}\), from m_parameters module)

• cohsesion_weak – upper cut-off value (\(\sigma_Y^{max}\), from m_parameters module)

• eps_strong – upper cut-off value (\(\sigma_Y^{max}\), from m_parameters module)

• eps_weak – upper cut-off value (\(\sigma_Y^{max}\), from m_parameters module)

Returns:

\[\begin{split}C = \begin{cases} C_0 &\text{if } \varepsilon_p < \varepsilon_0\\ C_0 + (C_\infty - C_0)\frac{\varepsilon_p - \varepsilon_0}{\varepsilon_\infty - \varepsilon_0} &\text{if } \varepsilon_0 < \varepsilon_p < \varepsilon_\infty\\ C_\infty & \text{if } \varepsilon_p \geq \varepsilon_\infty \end{cases}\end{split}\]

m_rheology.eta(Q, A, n, m, Temp, strain_rate_inv, stress_inv)

Evaluates the temperature- and stress-dependent viscosity for individual deformation mechanisms (dislocation creep, grain boundary slining and bassal slip).

Parameters:

• Q – activation energy for given mechanism

• A – exponential prefactor for given mechanism

• n – stress exponent for given mechanism

• m – grain size exponent for given mechanism

• Temp – temperature

• stress_invariant – second invariant of the deviatoric part of the Cauchy stress tensor

• strain_rate_invariant – second invariant of the strain rate tensor

Variables:

elasticity – Chooses between formulas depending whether elasticity is turned on or not.

Returns:

If elasticity = False, formula using the strain rate invariant is used

\[\eta_{i} = \frac{1}{3^{(n+1)/2n}\ 2^{(n-1)/n}} A^{-1/n}d^{m/n}\dot{\varepsilon}_{II}^{(1-n)/n}\exp\left( \frac{Q}{nRT} \right)\]

If elasticity = True, formula using the Cauchy stress invariant is used

\[\eta_{i} = \frac{1}{3^{(n+1)/2n}} A^{-1}d^{m}\dot{\sigma}_{II}^{1-n}\exp\left( \frac{Q}{RT} \right)\]

m_rheology.eta_BS(Temp, strain_rate_II, stress_II)

The same as m_rheology.eta_disl(), but for basal slip parameters.

m_rheology.eta_GBS(Temp, strain_rate_II, stress_II, eval_type)

The same as m_rheology.eta_disl(), but for grain boundary sliding parameters.

m_rheology.eta_diff(Temp)

Computes the viscosity resulting from the diffusion creep, taking into account the boundary and the volumetric diffusion \(\eta_{\rm diff}\).

Parameters:

Temp – temperature (\(T\))

Variables:

d_grain – grain size

Returns:

\[\eta_{diff} = \frac{1}{2 \cdot 3/2 \cdot 14 V_m} RTd^2 \left(D_v + \frac{\pi\delta D_b}{d} \right)^{-1}\]

where

\[\begin{split}D_{v/b} = D_{0,v/b}\exp\left(-\frac{Q}{RT} \right)\\\end{split}\]

m_rheology.eta_disl(Temp, strain_rate_II, stress_II, eval_type)

Computes the viscosity resulting from the dislocation creep \(\eta_{disl}\).

Parameters:

• Temp – temperature (\(T\))

• strain_rate_II – second invariant of the strain rate tensor (\(\dot{\varepsilon}_{II}\))

• stress_II – second invariant of the deviatoric part of the Cauchy stress tensor (\(\sigma_{II}\))

• eval_type – determines the method of the jump at the premelting temperature. String "mesh" returns conditional, "local" returns a simple if

Variables:

• n – stress coefficient

• p – grain size coefficient

• T_crit – temperature of premelting (\(T_{crit}\))

• Q_below – activation energy below the premelting temperature

• A_below – prefactor below the premelting temperature

• Q_above – activation energy above the premelting temperature

• A_above – prefactor above the premelting temperature

Returns:

calls the function m_rheology.eta() with the parameters of dislocation creep

m_rheology.eta_ductile(Temp, strain_rate, stress, xm, composition, step, Picard_iter, eval_type)

Evaluates the ductile viscosity based on the viscosity_type choice in m_parameters module.

Parameters:

• mesh – computational domain

• Temp – temperature

• xm – melt fraction

• stress_invariant – second invariant of the deviatoric part of the Cauchy stress tensor

• strain_rate_invariant – second invariant of the strain rate tensor

Variables:

viscosity_type – method for computing the viscosity, specified in the m_parameters module.

Returns:

• "constant" for constant viscosity of value \(\eta_0\).

\[\eta = \eta_0\]

• "temp-dep" for temperature-dependent viscosity

\[\eta(T) =\eta_0\cdot\exp\left[\frac{Q}{R}\left(\frac{1}{T} - \frac{1}{T_{\rm melt}}\right)\right]\]

• "GK_2001" for temperature and stress-dependent viscosity of ice-I

\[\eta(T, \sigma_{\rm II}) = \left(\frac{1}{\eta_{\rm diff}} + \frac{1}{\eta_{\rm disl}} + \frac{1}{\eta_{\rm BS} + \eta_{\rm GBS}} + \frac{1}{\eta_{\rm max} } \right)^{-1}\]

• "composition" for composition-dependent viscosity (used for selected benchmarks)

Tip

Custom viscosity function can be easily defined here.

m_rheology.eta_eff(p_k, Temp, strain_rate, stress, xm, composition, plastic_strain, step, Picard_iter, stress_dev_inv_k, shear_modulus, dt, sr_min)

Evaluates the ‘visco-plastic’ viscosity \(\eta_{vp}\).

Parameters:

• mesh – computational domain

• Temp – temperature

• xm – melt fraction

• stress_invariant – second invariant of the deviatoric part of the Cauchy stress tensor

• strain_rate_invariant – second invariant of the strain rate tensor

Variables:

eta_min_plast – lower cut-off value for plastic viscosity \(\eta_{p}\). Maintains reasonable viscosity contrast between tectonic faults and their surrounfings, see m_parameters.py.

Returns:

\[\eta_{vp} = \textrm{min}(\textrm{max}(\eta_p,\ \eta_p^{min}) ,\ \eta_v)\]

m_rheology.get_new_stress(mesh, Temp, xm, stress_dev_inv, stress_dev_inv_k, strain_rate_invariant, composition, shear_modulus, step, Picard_iter, dt)

Evaluates the visco-elastic stress assuming fully ductile viscosity.

Parameters:

• mesh – computational domain

• Temp – temperature

• xm – melt fraction

• stress_dev_inv – second invariant of the deviatoric part of the Cauchy stress tensor

• strain_rate_inv – second invariant of the strain rate tensor

• strain_rate_inv_k – second invariant of the strain rate tensor from the previous time step

• composition – material composition

• shear_modulus – shear modulus

• step – number of the time step (needed for 0th step for viscosity)

• Picard_iter – number of the Picard iteration (needed for 0th iteration for viscosity)

• dt – length of the time step

Returns:

If elasticity = False,

\[\sigma_{II}^{visc} = 2Z\eta_{visc}\dot{\varepsilon}_{II} + (1-Z)\sigma_{II}^k\]

If elasticity = False,

\[\sigma_{II}^{visc} = 2\eta_{visc}\dot{\varepsilon}_{II}\]

m_rheology.get_new_stress_iter(mesh, Temp, xm, stress_dev_inv, stress_dev_inv_k, strain_rate_inv, composition, shear_modulus, step, Picard_iter, dt)

Performs local visco-elasto-plastic iterations to ensure consistency between viscosity and stress during Picard iterations. Parameters are the same as in m_rheology.get_new_stress().

Caution

These iterations are needed every time when elasticity = True and the viscosity is stress-dependent.

m_rheology.max_function(a, b)

Numerical implementation of the maximum of two values. Extends the python max() to functions on the mesh.

Parameters:

• a – value or function on the mesh

• b – value or function on the mesh

Returns:

\[\text{max}(a,b) = \frac{a + b}{2} + \frac{\text{abs}(a + b)}{2}\]

m_rheology.plastic_strain_integration(mesh, tracers_in_cells, tracers, dt, yield_function, strain_rate)

Performs integration and healing of the plastic strain \(\varepsilon_p\) on markers.

Parameters:

• yield_function – yield function (\(F = \sigma_{II} - \sigma_Y\))

• dt – time step length (\(\Delta t\))

• strain_rate – second invariant of the strain rate tensor (\(\dot{\varepsilon}_{II}\))

• healing_timescale – healing time scale of the material (see healing_timescale in m_parameters.py)

Returns:

\[\begin{split}\varepsilon_p = \begin{cases} \dfrac{\varepsilon_p^k + \Delta t \dot{\varepsilon}_{II}}{1 + \Delta t / \tau} &\text{if } F = 0\\ \dfrac{\varepsilon_p^k}{1 + \Delta t / \tau} & \text{if } F < 0\\ \end{cases}\end{split}\]

m_rheology.sigma_yield(p_k, plastic_strain, composition)

Computes the yield stress \(\sigma_Y\).

Parameters:

• p_k – pressure (\(p\))

• plastic_strain – plastic strain (\(\varepsilon_p\))

• composition – material composition (\(C\))

• yield_stress_min – lower cut-off value (\(\sigma_Y^{min}\), from m_parameters module)

• yield_stress_max – upper cut-off value (\(\sigma_Y^{max}\), from m_parameters module)

Returns:

\[\sigma_Y = \textrm{min}(\textrm{max}(p\sin(\varphi) + C(\varepsilon_p)\cos(\varphi),\ \sigma_Y^{min}) ,\ \sigma_Y^{max})\]

m_rheology.strain_rate_II(v)

Computes the second invariant of the strain rate tensor.

Parameters:

v – velocity

Returns:

\[\dot{\varepsilon}_{II} = \sqrt{\frac{ \dot{\boldsymbol{\varepsilon}} : \dot{\boldsymbol{\varepsilon}}}{2}}\]

m_rheology.tensor_2nd_invariant(tensor)

Computes the second invariant of a tensor.

Parameters:

A – any tensor

Returns:

\[\dot{A}_{II} = \sqrt{\frac{ \dot{\boldsymbol{A}} : \dot{\boldsymbol{A}}}{2}}\]

m_rheology.tidal_heating(visc)

Computes the tidal heating rate \(H\).

Parameters:

visc – viscosity (\(\eta\))

Returns:

Defines \(\nu = \eta_{peak}/\eta\), where \(\eta_{peak}\) is the viscosity at which the maximum of the tidal heating occurs.

• If heating_model = "Maxwell"

  \[H = \frac{2 H_{max}}{\nu + 1/\nu}\]

• If heating_model = "Andrade"

  \[H = \frac{2 H_{max}}{\nu + 1/\nu}\]

m_rheology.z(visc, shear_modulus, dt)

Computes the viscoelastic parameter \(Z\).

Parameters:

• visc – viscosity (\(\eta\))

• dt – time step length (\(\Delta t\))

Returns:

\[Z = \frac{\Delta t}{\Delta t + \frac{\eta}{G}}\]