Approximations

approximations

This module provides classes that emulate physical approximations of fluid dynamics systems by exposing methods to calculate specific terms in the corresponding mathematical equations. Users instantiate the appropriate class by providing relevant parameters and pass the instance to other objects, such as solvers. Under the hood, G-ADOPT queries variables and methods from the approximation.

DeviatoricStressMixin

Provides the single definition of the deviatoric part of a gradient.

All logic about how the deviatoric part is taken lives in deviatoric_tensor_from_grad, keyed off the compressible property. stress_from_grad multiplies that tensor by the shear coefficient and adds whatever volumetric part the approximation carries. The full stress and the weak (SIPG) boundary terms in the momentum equation build on these two methods, so a subclass that changes the stress form (for example a true 2D compressible model, or one with a bulk modulus) overrides here.

We note that deviatoric_tensor_from_grad is purely kinematic and carries no material coefficient, so the same expression serves: 1 - the velocity gradient of a Stokes problem 2 - the incremental displacement gradient of a viscoelastic one 3 - the boundary jump tensor of the SIPG penalty.

Every caller supplies its own coefficient, for example self.mu, effective_viscosity(dt) in IncompressibleMaxwellApproximation.stress, \(\mu_0\) in MaxwellApproximation.deviatoric_stress.

deviatoric_tensor_from_grad(gradient)

Twice the deviatoric symmetric part of a gradient-like tensor.

For a gradient-like tensor \(G\) this returns

\[ 2\,\mathrm{sym}(G) \]

in the incompressible case, and

\[ 2\,\mathrm{sym}(G) - \tfrac{2}{3}\,\mathrm{tr}(G)\,I = 2\,\mathrm{dev}(\mathrm{sym}(G)) \]

in the compressible case, where \(\mathrm{dev}(A) = A - \tfrac{1}{3}\mathrm{tr}(A) I\). The two branches agree only where \(\mathrm{tr}(G) = 0\).

Note that the physical meaning of what is returned here depends on what the caller passes in. Cases are:

- With $G = \nabla u$ and $u$ a *velocity*, as in the Stokes problems,
  the result is twice the strain rate, with units of 1/s, and the
  caller multiplies it by a viscosity in Pa s.
- With $G = \nabla u$ and $u$ an *incremental displacement*, as in the
  viscoelastic and GIA approximations, the result is twice the strain
  and is dimensionless, and the caller multiplies it by a shear
  modulus in Pa.
- With $G = n \otimes w$, the outer product of the boundary normal
  with a jump in the weakly imposed solution, the result is the
  boundary jump tensor of the SIPG penalty and is neither a strain nor
  a strain rate.

The returned tensor therefore has no fixed units. Only the first two cases combine with a coefficient to give a stress, and only where that coefficient is \(\mu\) is the result \(\sigma(u)/\mu\).

Note: The factor \(\tfrac{2}{3}\) is fixed regardless of dimension, matching the stress in the geodynamics community. In 2D the deviator of a \(2 \times 2\) tensor subtracts \(\tfrac{1}{2}\mathrm{tr}(A) I\), so with the fixed \(\tfrac{2}{3}\) the compressible branch does not return a trace-free tensor in 2D.

Returns:

Type	Description
Expr	A UFL expression for twice the deviatoric symmetric part of
Expr	gradient.

Source code in g-adopt/gadopt/approximations.py

def deviatoric_tensor_from_grad(self, gradient: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
    r"""Twice the deviatoric symmetric part of a gradient-like tensor.

    For a gradient-like tensor $G$ this returns

    $$ 2\,\mathrm{sym}(G) $$

    in the incompressible case, and

    $$ 2\,\mathrm{sym}(G) - \tfrac{2}{3}\,\mathrm{tr}(G)\,I
       = 2\,\mathrm{dev}(\mathrm{sym}(G)) $$

    in the compressible case, where $\mathrm{dev}(A) = A -
    \tfrac{1}{3}\mathrm{tr}(A) I$. The two branches agree only where
    $\mathrm{tr}(G) = 0$.

    Note that the physical meaning of what is returned here depends on what the caller
    passes in. Cases are:

        - With $G = \nabla u$ and $u$ a *velocity*, as in the Stokes problems,
          the result is twice the strain rate, with units of 1/s, and the
          caller multiplies it by a viscosity in Pa s.
        - With $G = \nabla u$ and $u$ an *incremental displacement*, as in the
          viscoelastic and GIA approximations, the result is twice the strain
          and is dimensionless, and the caller multiplies it by a shear
          modulus in Pa.
        - With $G = n \otimes w$, the outer product of the boundary normal
          with a jump in the weakly imposed solution, the result is the
          boundary jump tensor of the SIPG penalty and is neither a strain nor
          a strain rate.

    The returned tensor therefore has no fixed units. Only the first two
    cases combine with a coefficient to give a stress, and only where that
    coefficient is $\mu$ is the result $\sigma(u)/\mu$.

    Note: The factor $\tfrac{2}{3}$ is fixed regardless of dimension,
    matching the stress in the geodynamics community. In 2D the deviator
    of a $2 \times 2$ tensor subtracts $\tfrac{1}{2}\mathrm{tr}(A) I$, so
    with the fixed $\tfrac{2}{3}$ the compressible branch does not return
    a trace-free tensor in 2D.

    Returns:
      A UFL expression for twice the deviatoric symmetric part of
      `gradient`.

    """
    # Twice the symmetric part of the gradient: $2 sym(G)$.
    deviatoric_tensor = 2 * sym(gradient)
    if self.compressible:
        # Remove the volumetric part, giving $2 dev(sym(G))$.
        deviatoric_tensor -= 2 / 3 * tr(gradient) * Identity(gradient.ufl_shape[0])
    return deviatoric_tensor

stress_from_grad(gradient)

The stress a gradient-like tensor produces, including any bulk part.

For a gradient-like tensor \(G\) this returns \(\mu\,A(G)\), with $A = $ deviatoric_tensor_from_grad. An approximation that also carries a volumetric response adds it by overriding this method.

There are two callers. With \(G = \nabla u\) this is the stress itself. With \(G = n \otimes w\), the outer product of the boundary normal with the jump in the weakly imposed velocity, it is the response the SIPG penalty term pairs against the test function, which is how the penalty picks up the same material coefficients as the stress it penalises.

Returns:

Type	Description
Expr	A UFL expression for the stress produced by gradient.

Source code in g-adopt/gadopt/approximations.py

def stress_from_grad(self, gradient: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
    r"""The stress a gradient-like tensor produces, including any bulk part.

    For a gradient-like tensor $G$ this returns $\mu\,A(G)$, with
    $A = $ `deviatoric_tensor_from_grad`. An approximation that also carries a
    volumetric response adds it by overriding this method.

    There are two callers. With $G = \nabla u$ this is the stress itself.
    With $G = n \otimes w$, the outer product of the boundary normal with
    the jump in the weakly imposed velocity, it is the response the SIPG
    penalty term pairs against the test function, which is how the penalty
    picks up the same material coefficients as the stress it penalises.

    Returns:
      A UFL expression for the stress produced by `gradient`.

    """
    return self.mu * self.deviatoric_tensor_from_grad(gradient)

BaseApproximation

Bases: DeviatoricStressMixin, ABC

Base class to provide expressions for the coupled Stokes and Energy system.

The basic assumption is that we are solving (to be extended when needed)

div(dev_stress) + grad p + buoyancy(T, p) * khat = 0
div(rho_continuity * u) = 0
rhocp DT/Dt + linearized_energy_sink(u) * T
  = div(kappa * grad(Tbar + T)) + energy_source(u)

where the following terms are provided by Approximation methods:

• linearized_energy_sink(u) = 0 (BA), Di * rhobar * alphabar * g * w (EBA), or Di * rhobar * alphabar * w (TALA/ALA)
• kappa() is diffusivity or conductivity depending on rhocp()
• Tbar is 0 or reference temperature profile (ALA)
• dev_stress depends on the compressible property (False or True):
  • if compressible then $\sigma = 2\mu\,\mathrm{sym}(\nabla u)
  • \tfrac{2}{3}\mu\,(\nabla\cdot u)\,I$
  • if not compressible then \(\sigma = 2\mu\,\mathrm{sym}(\nabla u)\) and rho_continuity is assumed to be 1

compressible abstractmethod property

Defines compressibility.

Returns:

Type	Description
bool	A boolean signalling if the governing equations are in compressible form.

Tbar abstractmethod property

Defines the reference temperature profile.

Returns:

Type	Description
Function	A Firedrake function for the reference temperature profile.

stress(u)

Defines the deviatoric stress \(\sigma(u)\).

Built from the single definition of the stress produced by a gradient-like tensor, so the incompressible form \(\sigma = 2\mu\,\mathrm{sym}(\nabla u)\) and the compressible form \(\sigma = 2\mu\,\mathrm{dev}(\mathrm{sym}(\nabla u))\) differ only through the compressible property.

Returns:

Type	Description
Expr	A UFL expression for the deviatoric stress.

Source code in g-adopt/gadopt/approximations.py

def stress(self, u: Function) -> ufl.core.expr.Expr:
    r"""Defines the deviatoric stress $\sigma(u)$.

    Built from the single definition of the stress produced by a
    gradient-like tensor, so the incompressible form
    $\sigma = 2\mu\,\mathrm{sym}(\nabla u)$ and the compressible form
    $\sigma = 2\mu\,\mathrm{dev}(\mathrm{sym}(\nabla u))$ differ only
    through the `compressible` property.

    Returns:
      A UFL expression for the deviatoric stress.

    """
    return self.stress_from_grad(grad(u))

buoyancy(p, T) abstractmethod

Defines the buoyancy force.

Returns:

Type	Description
Expr	A UFL expression for the buoyancy term (momentum source in gravity direction).

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def buoyancy(self, p: Function, T: Function) -> ufl.core.expr.Expr:
    """Defines the buoyancy force.

    Returns:
      A UFL expression for the buoyancy term (momentum source in gravity direction).

    """
    pass

rho_continuity() abstractmethod

Defines density.

Returns:

Type	Description
Expr	A UFL expression for density in the mass continuity equation.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def rho_continuity(self) -> ufl.core.expr.Expr:
    """Defines density.

    Returns:
      A UFL expression for density in the mass continuity equation.

    """
    pass

rhocp() abstractmethod

Defines the volumetric heat capacity.

Returns:

Type	Description
Expr	A UFL expression for the volumetric heat capacity in the energy equation.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def rhocp(self) -> ufl.core.expr.Expr:
    """Defines the volumetric heat capacity.

    Returns:
      A UFL expression for the volumetric heat capacity in the energy equation.

    """
    pass

kappa() abstractmethod

Defines thermal diffusivity.

Returns:

Type	Description
Expr	A UFL expression for thermal diffusivity.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def kappa(self) -> ufl.core.expr.Expr:
    """Defines thermal diffusivity.

    Returns:
      A UFL expression for thermal diffusivity.

    """
    pass

linearized_energy_sink(u) abstractmethod

Defines temperature-related sink terms.

Returns:

Type	Description
Expr	A UFL expression for temperature-related sink terms in the energy equation.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def linearized_energy_sink(self, u) -> ufl.core.expr.Expr:
    """Defines temperature-related sink terms.

    Returns:
      A UFL expression for temperature-related sink terms in the energy equation.

    """
    pass

energy_source(u) abstractmethod

Defines additional terms.

Returns:

Type	Description
Expr	A UFL expression for additional independent terms in the energy equation.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def energy_source(self, u) -> ufl.core.expr.Expr:
    """Defines additional terms.

    Returns:
      A UFL expression for additional independent terms in the energy equation.

    """
    pass

free_surface_terms(p, T, eta, theta_fs) abstractmethod

Defines free surface normal stress in the momentum equation and prefactor multiplying the free surface equation.

The normal stress depends on the density contrast across the free surface. Depending on the variable_rho_fs argument, this contrast either involves the interior reference density or the full interior density that accounts for changes in temperature and composition within the domain. For dimensional simulations, the user should specify delta_rho_fs, defined as the difference between the reference interior density and the exterior density. For non-dimensional simulations, the user should specify RaFS, the equivalent of the Rayleigh number for the density contrast across the free surface.

The prefactor multiplies the free surface equation to ensure the top-right and bottom-left corners of the block matrix remain symmetric. This is similar to rescaling eta -> eta_tilde in Kramer et al. (2012, see block matrix shown in Eq. 23).

Returns:

Type	Description
tuple[Expr]	A UFL expression for the free surface normal stress and a UFL expression for
tuple[Expr]	the free surface equation prefactor.

Source code in g-adopt/gadopt/approximations.py

@abc.abstractmethod
def free_surface_terms(self, p, T, eta, theta_fs) -> tuple[ufl.core.expr.Expr]:
    """Defines free surface normal stress in the momentum equation and prefactor
    multiplying the free surface equation.

    The normal stress depends on the density contrast across the free surface.
    Depending on the `variable_rho_fs` argument, this contrast either involves the
    interior reference density or the full interior density that accounts for
    changes in temperature and composition within the domain. For dimensional
    simulations, the user should specify `delta_rho_fs`, defined as the difference
    between the reference interior density and the exterior density. For
    non-dimensional simulations, the user should specify `RaFS`, the equivalent of
    the Rayleigh number for the density contrast across the free surface.

    The prefactor multiplies the free surface equation to ensure the top-right and
    bottom-left corners of the block matrix remain symmetric. This is similar to
    rescaling eta -> eta_tilde in Kramer et al. (2012, see block matrix shown in
    Eq. 23).

    Returns:
      A UFL expression for the free surface normal stress and a UFL expression for
      the free surface equation prefactor.

    """
    pass

BoussinesqApproximation(Ra, *, mu=1, rho=1, alpha=1, T0=0, g=1, RaB=0, delta_rho=1, kappa=1, H=0)

Bases: BaseApproximation

Expressions for the Boussinesq approximation.

Density variations are considered small and only affect the buoyancy term. Reference parameters are typically constant. Viscous dissipation is neglected (Di << 1).

Parameters:

Name	Type	Description	Default
Ra	Function | Number	Rayleigh number	required
mu	Function | Number	dynamic viscosity	1
rho	Function | Number	reference density	1
alpha	Function | Number	coefficient of thermal expansion	1
T0	Function | Number	reference temperature	0
g	Function | Number	gravitational acceleration	1
RaB	Function | Number	compositional Rayleigh number; product of the Rayleigh and buoyancy numbers	0
delta_rho	Function | Number	compositional density difference from the reference density	1
kappa	Function | Number	thermal diffusivity	1
H	Function | Number	internal heating rate	0

Note

The thermal diffusivity, gravitational acceleration, reference density, and coefficient of thermal expansion are normally kept at 1 when non-dimensionalised.

Source code in g-adopt/gadopt/approximations.py

def __init__(
    self,
    Ra: Function | Number,
    *,
    mu: Function | Number = 1,
    rho: Function | Number = 1,
    alpha: Function | Number = 1,
    T0: Function | Number = 0,
    g: Function | Number = 1,
    RaB: Function | Number = 0,
    delta_rho: Function | Number = 1,
    kappa: Function | Number = 1,
    H: Function | Number = 0,
):
    self.Ra = ensure_constant(Ra)
    self.mu = ensure_constant(mu)
    self.rho = ensure_constant(rho)
    self.alpha = ensure_constant(alpha)
    self.T0 = T0
    self.g = ensure_constant(g)
    self.kappa_ref = ensure_constant(kappa)
    self.RaB = RaB
    self.delta_rho = ensure_constant(delta_rho)
    self.H = ensure_constant(H)

ExtendedBoussinesqApproximation(Ra, Di, *, H=None, heating_weight=1, **kwargs)

Bases: BoussinesqApproximation

Expressions for the extended Boussinesq approximation.

Extends the Boussinesq approximation by including viscous dissipation and work against gravity (both scaled with Di).

Parameters:

Name	Type	Description	Default
Ra	Number	Rayleigh number	required
Di	Number	Dissipation number	required
mu		dynamic viscosity	required
H	Optional[Number]	volumetric heat production	None
heating_weight	Function | Number	scaling factor (Function or Number) applied to viscous dissipation and adiabatic heating terms; defaults to 1	1

Other Parameters:

Name	Type	Description
rho	Number	reference density
alpha	Number	coefficient of thermal expansion
T0	Function | Number	reference temperature
g	Number	gravitational acceleration
RaB	Number	compositional Rayleigh number; product of the Rayleigh and buoyancy numbers
delta_rho	Number	compositional density difference from the reference density
kappa	Number	thermal diffusivity

Note

The thermal diffusivity, gravitational acceleration, reference density, and coefficient of thermal expansion are normally kept at 1 when non-dimensionalised.

Source code in g-adopt/gadopt/approximations.py

def __init__(self, Ra: Number, Di: Number, *, H: Optional[Number] = None, heating_weight: Function | Number = 1, **kwargs):
    super().__init__(Ra, **kwargs)
    self.Di = Di
    self.H = H
    self.heating_weight = ensure_constant(heating_weight)

TruncatedAnelasticLiquidApproximation(Ra, Di, *, Tbar=0, cp=1, **kwargs)

Bases: ExtendedBoussinesqApproximation

Truncated Anelastic Liquid Approximation

Compressible approximation. Excludes linear dependence of density on pressure.

Parameters:

Name	Type	Description	Default
Ra	Number	Rayleigh number	required
Di	Number	Dissipation number	required
Tbar	Function | Number	reference temperature. In the diffusion term we use Tbar + T (i.e. T is the pertubartion)	0
cp	Function | Number	reference specific heat at constant pressure	1

Other Parameters:

Name	Type	Description
rho	Number	reference density
alpha	Number	reference thermal expansion coefficient
T0	Function | Number	reference temperature
g	Number	gravitational acceleration
RaB	Number	compositional Rayleigh number; product of the Rayleigh and buoyancy numbers
delta_rho	Number	compositional density difference from the reference density
kappa	Number	diffusivity
mu	Number	viscosity used in viscous dissipation
H	Number	volumetric heat production

Note

Other keyword arguments may be depth-dependent, but default to 1 if not supplied.

Source code in g-adopt/gadopt/approximations.py

def __init__(self,
             Ra: Number,
             Di: Number,
             *,
             Tbar: Function | Number = 0,
             cp: Function | Number = 1,
             **kwargs):
    super().__init__(Ra, Di, **kwargs)
    self.Tbar = Tbar
    self.cp = cp

AnelasticLiquidApproximation(Ra, Di, *, chi=1, gamma0=1, cp0=1, cv0=1, **kwargs)

Bases: TruncatedAnelasticLiquidApproximation

Anelastic Liquid Approximation

Compressible approximation. Includes linear dependence of density on pressure.

Parameters:

Name	Type	Description	Default
Ra	Number	Rayleigh number	required
Di	Number	Dissipation number	required
chi	Function | Number	reference isothermal compressibility	1
gamma0	Function | Number	Gruneisen number (in pressure-dependent buoyancy term)	1
cp0	Function | Number	specific heat at constant pressure, reference for entire Mantle (in pressure-dependent buoyancy term)	1
cv0	Function | Number	specific heat at constant volume, reference for entire Mantle (in pressure-dependent buoyancy term)	1

Other Parameters:

Name	Type	Description
rho	Number	reference density
alpha	Number	reference thermal expansion coefficient
T0	Function | Number	reference temperature
g	Number	gravitational acceleration
RaB	Number	compositional Rayleigh number; product of the Rayleigh and buoyancy numbers
delta_rho	Number	compositional density difference from the reference density
kappa	Number	diffusivity
mu	Number	viscosity used in viscous dissipation
H	Number	volumetric heat production
Tbar	Number	reference temperature. In the diffusion term we use Tbar + T (i.e. T is the pertubartion)
cp	Number	reference specific heat at constant pressure

Source code in g-adopt/gadopt/approximations.py

def __init__(self,
             Ra: Number,
             Di: Number,
             *,
             chi: Function | Number = 1,
             gamma0: Function | Number = 1,
             cp0: Function | Number = 1,
             cv0: Function | Number = 1,
             **kwargs):
    super().__init__(Ra, Di, **kwargs)
    # Dynamic pressure contribution towards buoyancy
    self.chi = chi
    self.gamma0, self.cp0, self.cv0 = gamma0, cp0, cv0

BaseGIAApproximation(density, shear_modulus, viscosity, *, g=1, B_mu=1)

Bases: DeviatoricStressMixin

Base class for viscoelasticity assuming small displacements.

By assuming a small displacement with respect to mantle thickness we can linearise the problem, introducing a perturbation away from a reference state.

For background and derivation of the formulation please see the equations and references provided in Scott et al 2025.

Automated forward and adjoint modelling of viscoelastic deformation of the solid Earth. Scott, W.; Hoggard, M.; Duvernay, T.; Ghelichkhan, S.; Gibson, A.; Roberts, D.; Kramer, S. C.; and Davies, D. R. EGUsphere, 2025: 1–43. 2025.

Parameters:

Name	Type	Description	Default
density	Function | Number	density of the reference state - assumed to be hydrostatic	required
shear_modulus	Function | Number | list	shear modulus	required
viscosity	Function | Number | list	viscosity	required
g	Function | Number	gravitational acceleration	1
B_mu	Function | Number	Nondimensional number describing ratio of buoyancy to elastic shear strength used for nondimensionalisation. $$ B_{\mu} = rac{ar{ ho} ar{g} L}{ar{\mu}}$, where \(ar{ ho}\) is a characteristic density scale (kg / m^3), \(ar{g}\) is a characteristic gravity scale (m / s^2), \(L\) is a characteristic length scale, often Mantle depth (m), \(\mu\) is a characteristic shear modulus (Pa).	1

Source code in g-adopt/gadopt/approximations.py

def __init__(
    self,
    density: Function | Number,
    shear_modulus: Function | Number | list,
    viscosity: Function | Number | list,
    *,
    g: Function | Number = 1,
    B_mu: Function | Number = 1,
):
    self.density = ensure_constant(density)
    self.shear_modulus = ensure_constant(shear_modulus)
    self.viscosity = ensure_constant(viscosity)
    self.g = ensure_constant(g)
    self.B_mu = ensure_constant(B_mu)

IncompressibleMaxwellApproximation(density, shear_modulus, viscosity, **kwargs)

Bases: BaseGIAApproximation

Incompressible Maxwell rheology via the incremental displacement formulation.

This class implements incompressible Maxwell rheology similar to the approach by Zhong et al. 2003. The linearised problem is cast in terms of incremental displacement, i.e. velocity * dt where dt is the timestep. This produces a mixed Stokes system for incremental displacement and pressure which can be solved in the same way as mantle convection (where unknowns are velocity and pressure), with a modfied viscosity and stress term accounting for the deviatoric stress at the previous timestep.

Zhong, Shijie, Archie Paulson, and John Wahr. "Three-dimensional finite-element modelling of Earth’s viscoelastic deformation: effects of lateral variations lithospheric thickness." Geophysical Journal International 155.2 (2003): 679-695.

N.b. that the implementation currently assumes all terms are dimensional.

Parameters:

Name	Type	Description	Default
density	Function | Number	background density	required
shear_modulus	Function | Number	shear modulus	required
viscosity	Function | Number	viscosity	required
g		gravitational acceleration	required

Source code in g-adopt/gadopt/approximations.py

def __init__(
    self,
    density: Function | Number,
    shear_modulus: Function | Number,
    viscosity: Function | Number,
    **kwargs,
):

    warn(
        '''The IncompressibleMaxwellApproximation is being phased out of G-ADOPT.
        We recommend using `MaxwellApproximation` together with the
        `InternalVariableSolver` for viscoelastic applications in G-ADOPT.
    ''')

    super().__init__(density, shear_modulus, viscosity, **kwargs)
    self.maxwell_time = self.viscosity / self.shear_modulus

InternalVariableApproximation(bulk_modulus, density, shear_modulus, viscosity, *, bulk_shear_ratio=1, exponent=1, transition_stress=1, background_stress=0, **kwargs)

Bases: BaseGIAApproximation

Base class for viscoelastic stress relaxation via internal variables.

This class implements compressible viscoelasticity following the formulation adopted by Al-Attar and Tromp (2014) and Crawford et al. (2017, 2018), in which viscoelastic constitutive equations are expressed in integral form and reformulated using so-called internal variables. Conceptually, this approach consists of a set of elements with different shear relaxation timescales, arranged in parallel. This formulation provides a compact, flexible and convenient means to incorporate transient rheology into viscoelastic deformation models: using a single internal variable is equivalent to a simple Maxwell material; two correspond to a Burgers model with two characteristic relaxation frequencies; and using a series of internal variables permits approximation of a continuous range of relaxation timescales for more complicated rheologies.

For more information regarding the specific implementation in G-ADOPT please see Scott et al. 2025.

Al-Attar, David, and Jeroen Tromp. "Sensitivity kernels for viscoelastic loading based on adjoint methods." Geophysical Journal International 196.1 (2014): 34-77.

Crawford, O., Al-Attar, D., Tromp, J., & Mitrovica, J. X. (2016). Forward and inverse modelling of post-seismic deformation. Geophysical Journal International, ggw414.

Crawford, O., Al-Attar, D., Tromp, J., Mitrovica, J. X., Austermann, J., & Lau, H. C. (2018). Quantifying the sensitivity of post-glacial sea level change to laterally varying viscosity. Geophysical journal international, 214(2), 1324-1363.

Automated forward and adjoint modelling of viscoelastic deformation of the solid Earth. Scott, W.; Hoggard, M.; Duvernay, T.; Ghelichkhan, S.; Gibson, A.; Roberts, D.; Kramer, S. C.; and Davies, D. R. EGUsphere, 2025: 1–43. 2025.

Parameters:

Name	Type	Description	Default
bulk_modulus	Function | Number	bulk modulus	required
density	Function | Number	density of the reference state - assumed to be hydrostatic	required
shear_modulus	Function | Number | list	shear modulus	required
viscosity	Function | Number | list	viscosity	required
bulk_shear_ratio	Function | Number	Ratio of bulk to shear modulus	1
exponent	Function | Number	Power-law stress exponent n. Defaults to 1, which recovers Newtonian (linear) viscosity with no special casing. Set n > 1 for composite (diffusion + dislocation) creep. The power-law factor is always evaluated; for n=1 it is identically 1.	1
transition_stress	Function | Number	Stress level (in sqrt(J2) units, consistent with second_stress_invariant) at which diffusion and dislocation creep contribute equally. Defaults to 1; the value is irrelevant when exponent=1.	1
background_stress	Function | Number	Background stress magnitude (sqrt(J2)) representing a pre-existing stress state, e.g. from mantle convection. Defaults to 0. The total stress entering the power-law factor is sigma_GIA + background_stress, a scalar combination that assumes the two fields are approximately aligned.	0
g		gravitational acceleration	required
B_mu		Nondimensional number describing ratio of buoyancy to elastic shear strength used for nondimensionalisation. $$ B_{\mu} = rac{ar{ ho} ar{g} L}{ar{\mu}}$, where \(ar{ ho}\) is a characteristic density scale (kg / m^3), \(ar{g}\) is a characteristic gravity scale (m / s^2), \(L\) is a characteristic length scale, often Mantle depth (m), \(\mu\) is a characteristic shear modulus (Pa).	required

Source code in g-adopt/gadopt/approximations.py

def __init__(
    self,
    bulk_modulus: Function | Number,
    density: Function | Number,
    shear_modulus: Function | Number | list,
    viscosity: Function | Number | list,
    *,
    bulk_shear_ratio: Function | Number = 1,
    exponent: Function | Number = 1,
    transition_stress: Function | Number = 1,
    background_stress: Function | Number = 0,
    **kwargs,
):
    self.bulk_modulus = ensure_constant(bulk_modulus)
    self.bulk_shear_ratio = ensure_constant(bulk_shear_ratio)

    # Convert viscosity and shear modulus to lists in the case
    # for Maxwell Rheology where there is only one internal variable
    # and hence only one pair of viscosity and shear modulus fields
    if not isinstance(shear_modulus, list):
        shear_modulus = [shear_modulus]
    if not isinstance(viscosity, list):
        viscosity = [viscosity]

    if len(viscosity) != len(shear_modulus):
        raise ValueError(
            "Length of viscosity and shear modulus lists must be consistent"
        )

    super().__init__(density, shear_modulus, viscosity, **kwargs)
    self.maxwell_times = [
        ensure_constant(visc / mu)
        for visc, mu in zip(self.viscosity, self.shear_modulus)
    ]
    # GIA literature uses `mu` for the shear modulus. Here, `mu` is the
    # solver-selected shear scale for the Nitsche penalty and solver context.
    # The material lists contain one entry per Maxwell element, and `mu0` is
    # their total elastic shear modulus.
    self.mu0 = ensure_constant(sum(self.shear_modulus))

    # Power law arguments
    self.exponent = exponent
    self.transition_stress = transition_stress
    self.background_stress = background_stress

deviatoric_strain(u)

The deviatoric strain \(\mathrm{dev}(\mathrm{sym}(\nabla u))\).

This is half of deviatoric_tensor_from_grad, which for this compressible approximation is \(2\,\mathrm{sym}(\nabla u) - \tfrac{2}{3}(\nabla \cdot u) I\). Sharing that one definition keeps the strain the internal variables relax towards and the stress the momentum equation assembles from the same expression. Here \(u\) is an incremental displacement, so the result is a strain and is dimensionless.

N.b. the fixed \(\tfrac{2}{3}\) makes this slightly inconsistent in 2D, for the reason given in DeviatoricStressMixin. The analytical tests in tests/viscoelastic_internal_variable/ are set up assuming 3D geometry.

Returns:

Type	Description
Expr	A UFL expression for the deviatoric strain.

Source code in g-adopt/gadopt/approximations.py

def deviatoric_strain(self, u: Function) -> ufl.core.expr.Expr:
    r"""The deviatoric strain $\mathrm{dev}(\mathrm{sym}(\nabla u))$.

    This is half of `deviatoric_tensor_from_grad`, which for this
    compressible approximation is $2\,\mathrm{sym}(\nabla u)
    - \tfrac{2}{3}(\nabla \cdot u) I$. Sharing that one definition keeps
    the strain the internal variables relax towards and the stress the
    momentum equation assembles from the same expression. Here $u$ is an
    incremental displacement, so the result is a strain and is
    dimensionless.

    N.b. the fixed $\tfrac{2}{3}$ makes this slightly inconsistent in 2D,
    for the reason given in `DeviatoricStressMixin`. The analytical tests
    in tests/viscoelastic_internal_variable/ are set up assuming 3D
    geometry.

    Returns:
      A UFL expression for the deviatoric strain.

    """
    return 0.5 * self.deviatoric_tensor_from_grad(grad(u))

effective_viscosity(dt)

Shear coefficient of the time-step problem after substitution.

The backward-Euler substitution of the internal variables gives a displacement problem with shear coefficient

\[ \eta_{eff} = \sum_i \frac{\eta_i}{\tau_i + \Delta t} = \sum_i \frac{\mu_i}{1 + \Delta t/\tau_i}, \]

a modulus despite the name. InternalVariableSolver uses this value for its substituted displacement operator. CoupledInternalVariableSolver uses mu0 for its elastic displacement block.

Parameters:

Name	Type	Description	Default
dt	float	the time step.	required

Returns:

Type	Description
Expr	A UFL expression for the effective shear coefficient.

Source code in g-adopt/gadopt/approximations.py

def effective_viscosity(self, dt: float) -> ufl.core.expr.Expr:
    r"""Shear coefficient of the time-step problem after substitution.

    The backward-Euler substitution of the internal variables gives a
    displacement problem with shear coefficient

    $$ \eta_{eff} = \sum_i \frac{\eta_i}{\tau_i + \Delta t}
                 = \sum_i \frac{\mu_i}{1 + \Delta t/\tau_i}, $$

    a modulus despite the name. `InternalVariableSolver` uses this value for
    its substituted displacement operator. `CoupledInternalVariableSolver`
    uses `mu0` for its elastic displacement block.

    Arguments:
      dt: the time step.

    Returns:
      A UFL expression for the effective shear coefficient.

    """

    eta_eff = 0
    for eta, maxwell_time in zip(self.viscosity, self.maxwell_times):
        eta_eff += eta / (maxwell_time + dt)
    return eta_eff

stress_from_grad(gradient)

The stress a gradient-like tensor produces, deviatoric plus bulk.

This approximation is compressible with a finite bulk modulus, so the response to a gradient-like tensor \(G\) has a volumetric part on top of the deviatoric one,

\[ \mu\,A(G) + \kappa_r\,\kappa\,\mathrm{tr}(G)\,I, \]

with \(\kappa\) the bulk modulus and \(\kappa_r\) the bulk-to-shear ratio. With \(G = \nabla u\) the volumetric part is \(\kappa_r \kappa (\nabla \cdot u) I\). With \(G = n \otimes w\) it is the volumetric part of the SIPG penalty, which the boundary jump in the normal displacement must feel for the same reason the interior divergence does.

The solver assigns \(\mu\) before assembly. InternalVariableSolver uses the effective shear coefficient of its substituted problem. CoupledInternalVariableSolver uses \(\mu_0\) to match its elastic displacement block.

Returns:

Type	Description
Expr	A UFL expression for the stress produced by gradient.

Source code in g-adopt/gadopt/approximations.py

def stress_from_grad(self, gradient: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
    r"""The stress a gradient-like tensor produces, deviatoric plus bulk.

    This approximation is compressible with a finite bulk modulus, so the
    response to a gradient-like tensor $G$ has a volumetric part on top of
    the deviatoric one,

    $$ \mu\,A(G) + \kappa_r\,\kappa\,\mathrm{tr}(G)\,I, $$

    with $\kappa$ the bulk modulus and $\kappa_r$ the bulk-to-shear ratio.
    With $G = \nabla u$ the volumetric part is
    $\kappa_r \kappa (\nabla \cdot u) I$. With $G = n \otimes w$ it is the
    volumetric part of the SIPG penalty, which the boundary jump in the
    normal displacement must feel for the same reason the interior
    divergence does.

    The solver assigns $\mu$ before assembly. `InternalVariableSolver` uses
    the effective shear coefficient of its substituted problem.
    `CoupledInternalVariableSolver` uses $\mu_0$ to match its elastic
    displacement block.

    Returns:
      A UFL expression for the stress produced by `gradient`.

    """
    volumetric = self.bulk_shear_ratio * self.bulk_modulus * tr(gradient)
    return (
        super().stress_from_grad(gradient)
        + volumetric * Identity(gradient.ufl_shape[0])
    )

stress(u, **kwargs)

The full stress, volumetric part plus relaxing deviatoric part.

\[ \sigma = \kappa_r\,\kappa\,(\nabla \cdot u)\,I + 2\mu_0\,\mathrm{dev}\,\varepsilon(u) - \sum_i 2\mu_i m_i. \]

The shear coefficient of the deviatoric part is the elastic modulus \(\mu_0 = \sum_i \mu_i\); the internal variables \(m_i\) carry the relaxation away from that elastic response. This is deliberately not written through stress_from_grad, whose shear coefficient is the solver-selected scale mu, not always \(\mu_0\).

Returns:

Type	Description
Expr	A UFL expression for the full stress.

Source code in g-adopt/gadopt/approximations.py

def stress(self, u, **kwargs) -> ufl.core.expr.Expr:
    r"""The full stress, volumetric part plus relaxing deviatoric part.

    $$ \sigma = \kappa_r\,\kappa\,(\nabla \cdot u)\,I
       + 2\mu_0\,\mathrm{dev}\,\varepsilon(u)
       - \sum_i 2\mu_i m_i. $$

    The shear coefficient of the deviatoric part is the elastic modulus
    $\mu_0 = \sum_i \mu_i$; the internal variables $m_i$ carry the
    relaxation away from that elastic response. This is deliberately not
    written through `stress_from_grad`, whose shear coefficient is the
    solver-selected scale `mu`, not always $\mu_0$.

    Returns:
      A UFL expression for the full stress.

    """
    internal_variables = kwargs.get("internal_variables", None)
    if internal_variables is None:
        raise KeyError(
            f"The internal_variables kwarg must be provided for stress() in {self.__class__.__name__}"
        )

    div_u = div(u) * Identity(len(u))
    stress = self.bulk_shear_ratio * self.bulk_modulus * div_u
    stress += self.deviatoric_stress(u, internal_variables)
    return stress

deviatoric_stress(u, internal_variables)

The deviatoric stress, elastic response minus the relaxed part.

\[ 2\mu_0\,\mathrm{dev}\,\varepsilon(u) - \sum_i 2\mu_i m_i \]

Built on deviatoric_tensor_from_grad, which supplies \(2\,\mathrm{dev}\,\varepsilon(u)\), so the elastic part of the stress and the strain the internal variables relax towards share one definition.

Returns:

Type	Description
	A UFL expression for the deviatoric stress.

Source code in g-adopt/gadopt/approximations.py

def deviatoric_stress(self, u, internal_variables):
    r"""The deviatoric stress, elastic response minus the relaxed part.

    $$ 2\mu_0\,\mathrm{dev}\,\varepsilon(u) - \sum_i 2\mu_i m_i $$

    Built on `deviatoric_tensor_from_grad`, which supplies
    $2\,\mathrm{dev}\,\varepsilon(u)$, so the elastic part of the stress
    and the strain the internal variables relax towards share one
    definition.

    Returns:
      A UFL expression for the deviatoric stress.

    """
    dev_stress = self.mu0 * self.deviatoric_tensor_from_grad(grad(u))
    for mu, m in zip(self.shear_modulus, internal_variables):
        dev_stress -= 2 * mu * m
    return dev_stress

second_stress_invariant(dev_stress)

Return the second deviatoric stress invariant sqrt(J2).

Computes sqrt(J2) = sqrt(0.5 * s_ij * s_ij), the standard second invariant used in mantle rheology power-law formulations. The transition_stress parameter throughout this class must be supplied in the same units, i.e. as a sqrt(J2) value.

Source code in g-adopt/gadopt/approximations.py

def second_stress_invariant(self, dev_stress):
    """Return the second deviatoric stress invariant sqrt(J2).

    Computes sqrt(J2) = sqrt(0.5 * s_ij * s_ij), the standard second
    invariant used in mantle rheology power-law formulations. The
    `transition_stress` parameter throughout this class must be supplied
    in the same units, i.e. as a sqrt(J2) value.
    """
    return sqrt(0.5 * inner(dev_stress, dev_stress) + 1e-16)

power_law_factor(dev_stress)

Return the composite creep viscosity correction factor.

For composite (diffusion + dislocation) creep the effective viscosity is

eta_eff(sigma) = eta_diff / (1 + (sigma / sigma*)^(n-1))

where sigma* is the transition stress (in sqrt(J2) units, consistent with second_stress_invariant) at which both mechanisms contribute equally. The input viscosity field is assumed to represent eta_eff evaluated at background_stress (e.g. from a mantle convection simulation). The factor returned here corrects that value to the total stress sigma_GIA + sigma_bg:

f = (1 + a^(n-1)) / (1 + b^(n-1))

where a = sigma_bg / sigma and b = (sigma_GIA + sigma_bg) / sigma.

Note: sigma_GIA and sigma_bg are combined as scalar magnitudes. This assumes the two stress fields are approximately aligned; it is an accepted simplification in GIA power-law literature but should be borne in mind when interpreting results.

Note on UFL evaluation: UFL evaluates 0^0 = 0 rather than the mathematically correct value of 1. A conditional is used to guard the exponentiation when the base is zero so that n=1 (Newtonian) layers correctly return f=1 regardless of background_stress.

Source code in g-adopt/gadopt/approximations.py

def power_law_factor(self, dev_stress):
    """Return the composite creep viscosity correction factor.

    For composite (diffusion + dislocation) creep the effective viscosity is

        eta_eff(sigma) = eta_diff / (1 + (sigma / sigma*)^(n-1))

    where sigma* is the transition stress (in sqrt(J2) units, consistent with
    second_stress_invariant) at which both mechanisms contribute equally.
    The input viscosity field is assumed to represent eta_eff evaluated at
    `background_stress` (e.g. from a mantle convection simulation). The factor
    returned here corrects that value to the total stress sigma_GIA + sigma_bg:

        f = (1 + a^(n-1)) / (1 + b^(n-1))

    where a = sigma_bg / sigma* and b = (sigma_GIA + sigma_bg) / sigma*.

    Note: sigma_GIA and sigma_bg are combined as scalar magnitudes. This
    assumes the two stress fields are approximately aligned; it is an accepted
    simplification in GIA power-law literature but should be borne in mind
    when interpreting results.

    Note on UFL evaluation: UFL evaluates 0^0 = 0 rather than the
    mathematically correct value of 1. A conditional is used to guard the
    exponentiation when the base is zero so that n=1 (Newtonian) layers
    correctly return f=1 regardless of background_stress.
    """
    sigma = self.second_stress_invariant(dev_stress)
    n = self.exponent
    a = self.background_stress / self.transition_stress
    b = (sigma + self.background_stress) / self.transition_stress
    # Guard against 0^0: UFL evaluates 0^0 as 0 rather than 1. This
    # only arises when n=1 (exponent n-1 = 0); for n>1, a zero base
    # raised to a positive power correctly evaluates to 0. Conditioning
    # on the exponent rather than the base avoids incorrectly setting
    # a^(n-1) = 1 for n>1 when background_stress = 0.
    a_pow = ufl.conditional(ufl.eq(n - 1, 0), ufl.as_ufl(1), a ** (n - 1))
    b_pow = ufl.conditional(ufl.eq(n - 1, 0), ufl.as_ufl(1), b ** (n - 1))
    return (1 + a_pow) / (1 + b_pow)

QuasiCompressibleInternalVariableApproximation(bulk_modulus, density, shear_modulus, viscosity, *, bulk_shear_ratio=1, **kwargs)

Bases: InternalVariableApproximation

Quasi compressible viscoelasticty via internal variables

N.b. this class neglects two key terms: compressible buoyancy and the volume integral after integrating the advection of hydrostatic prestress term by parts. This allows us to reproduce simplified analytical cases such as the Cathles 2024 benchmark in g-adopt/test/viscoelastic_internal_variable/ where these terms are not accounted for in the analytical solution.

For more information on the formulation see the documentation of the parent class InternalVariableApproximation

Parameters:

Name	Type	Description	Default
bulk_modulus	Function | Number	bulk modulus	required
density	Function | Number	density of the reference state - assumed to be hydrostatic	required
shear_modulus	Function | Number | list[Function | Number]	shear modulus	required
viscosity	Function | Number | list[Function | Number]	viscosity	required
bulk_shear_ratio	Function | Number	Ratio of bulk to shear modulus	1
exponent		Power-law stress exponent n (default 1 = Newtonian)	required
transition_stress		Transition stress in sqrt(J2) units (default 1, irrelevant for n=1)	required
background_stress		Background stress magnitude in sqrt(J2) units (default 0)	required
g		gravitational acceleration	required
B_mu		Nondimensional number describing ratio of buoyancy to elastic shear strength used for nondimensionalisation. $$ B_{\mu} = rac{ar{ ho} ar{g} L}{ar{\mu}}$, where \(ar{ ho}\) is a characteristic density scale (kg / m^3), \(ar{g}\) is a characteristic gravity scale (m / s^2), \(L\) is a characteristic length scale, often Mantle depth (m), \(\mu\) is a characteristic shear modulus (Pa).	required

Source code in g-adopt/gadopt/approximations.py

def __init__(
    self,
    bulk_modulus: Function | Number,
    density: Function | Number,
    shear_modulus: Function | Number | list[Function | Number],
    viscosity: Function | Number | list[Function | Number],
    *,
    bulk_shear_ratio: Function | Number = 1,
    **kwargs,
):

    warn(
        '''The QuasiCompressibleInternalVariableApproximation has specifically
        been designed to reproduce the simplified analytical benchmarks of
        Cathles 2024, e.g. see `g-adopt/tests/viscoelastic_internal_variable`.
        We recommend using `MaxwellApproximation` together with the
        `InternalVariableSolver` for viscoelastic applications in G-ADOPT.
    ''')

    super().__init__(
        bulk_modulus,
        density,
        shear_modulus,
        viscosity,
        bulk_shear_ratio=bulk_shear_ratio,
        **kwargs,
    )

CompressibleInternalVariableApproximation(bulk_modulus, density, shear_modulus, viscosity, *, bulk_shear_ratio=1, **kwargs)

Bases: InternalVariableApproximation

Compressible viscoelastic rheology via the internal variable formulation.

For more information on the formulation see the documentation of the parent class InternalVariableApproximation.

Parameters:

Name	Type	Description	Default
bulk_modulus	Function | Number	bulk modulus	required
density	Function | Number	density of the reference state - assumed to be hydrostatic	required
shear_modulus	Function | Number | list[Function | Number]	shear modulus	required
viscosity	Function | Number | list[Function | Number]	viscosity	required
bulk_shear_ratio	Function | Number	Ratio of bulk to shear modulus	1
exponent		Power-law stress exponent n (default 1 = Newtonian)	required
transition_stress		Transition stress in sqrt(J2) units (default 1, irrelevant for n=1)	required
background_stress		Background stress magnitude in sqrt(J2) units (default 0)	required
g		gravitational acceleration	required
B_mu		Nondimensional number describing ratio of buoyancy to elastic shear strength used for nondimensionalisation. $$ B_{\mu} = rac{ar{ ho} ar{g} L}{ar{\mu}}$, where \(ar{ ho}\) is a characteristic density scale (kg / m^3), \(ar{g}\) is a characteristic gravity scale (m / s^2), \(L\) is a characteristic length scale, often Mantle depth (m), \(\mu\) is a characteristic shear modulus (Pa).	required

Source code in g-adopt/gadopt/approximations.py

def __init__(
    self,
    bulk_modulus: Function | Number,
    density: Function | Number,
    shear_modulus: Function | Number | list[Function | Number],
    viscosity: Function | Number | list[Function | Number],
    *,
    bulk_shear_ratio: Function | Number = 1,
    **kwargs,
):

    super().__init__(
        bulk_modulus,
        density,
        shear_modulus,
        viscosity,
        bulk_shear_ratio=bulk_shear_ratio,
        **kwargs,
    )

MaxwellApproximation(bulk_modulus, density, shear_modulus, viscosity, **kwargs)

Bases: CompressibleInternalVariableApproximation

Maxwell viscoelastic rheology.

This is a helper class to simplify setting up (compressible) Maxwell rheology using the internal variable approach. For more references on the method please refer to the documentation of the parent class.

Parameters:

Name	Type	Description	Default
bulk_modulus	Function | Number	bulk modulus	required
density	Function | Number	background density	required
shear_modulus	Function | Number	shear modulus	required
viscosity	Function | Number	viscosity	required
bulk_shear_ratio		Ratio of bulk to shear modulus	required
g		gravitational acceleration	required
B_mu		Nondimensional number describing ratio of buoyancy to elastic shear strength used for nondimensionalisation. $ B_{\mu} = rac{ar{ ho} ar{g} L}{ar{\mu}}$, where \(ar{ ho}\) is a characteristic density scale (kg / m^3), \(ar{g}\) is a characteristic gravity scale (m / s^2), \(L\) is a characteristic length scale, often Mantle depth (m), \(\mu\) is a characteristic shear modulus (Pa).	required

Source code in g-adopt/gadopt/approximations.py

def __init__(
    self,
    bulk_modulus: Function | Number,
    density: Function | Number,
    shear_modulus: Function | Number,
    viscosity: Function | Number,
    **kwargs,
):
    super().__init__(bulk_modulus, density, shear_modulus, viscosity, **kwargs)