Skip to content

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.

BaseApproximation

Bases: 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 dev_stress = mu * [sym(grad(u) - 2/3 div(u)]
    • if not compressible then dev_stress = mu * sym(grad(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) abstractmethod

Defines the deviatoric stress.

Returns:

Type Description
Expr

A UFL expression for the deviatoric stress.
Source code in g-adopt/gadopt/approximations.py 
64
65
66
67
68
69
70
71
72
@abc.abstractmethod
def stress(self, u: Function) -> ufl.core.expr.Expr:
    """Defines the deviatoric stress.

    Returns:
      A UFL expression for the deviatoric stress.

    """
    pass

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 
74
75
76
77
78
79
80
81
82
@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 
84
85
86
87
88
89
90
91
92
@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 
 94
 95
 96
 97
 98
 99
100
101
102
@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 
104
105
106
107
108
109
110
111
112
@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 
125
126
127
128
129
130
131
132
133
@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 
135
136
137
138
139
140
141
142
143
@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 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
@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 
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
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, **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

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 
292
293
294
295
def __init__(self, Ra: Number, Di: Number, *, H: Optional[Number] = None, **kwargs):
    super().__init__(Ra, **kwargs)
    self.Di = Di
    self.H = H

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 
346
347
348
349
350
351
352
353
354
355
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 
400
401
402
403
404
405
406
407
408
409
410
411
412
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)

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 
450
451
452
453
454
455
456
457
458
459
460
461
462
463
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 
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
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

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

Bases: BaseGIAApproximation

Quasi compressible viscoelasticty 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.

This class implements the substitution method where the time-dependent internal variable equation is substituted into the momentum equation assuming a Backward Euler time discretisation. Therefore, the displacement field is the only unknown. For more information regarding the specific implementation in G-ADOPT please see Scott et al. 2025.

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 /tests/glacial_isostatic_adjustment/iv_ve_fs.py

where we need to remove these effect.

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.

Source code in g-adopt/gadopt/approximations.py 
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
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,
    **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)
    ]
    self.mu0 = ensure_constant(sum(self.shear_modulus))

effective_viscosity(dt)

Effective viscosity used to impose boundary conditions on displacement weakly through a Nitsche penalty term in the viscosity_term of momentum_equation.py

Source code in g-adopt/gadopt/approximations.py 
657
658
659
660
661
662
663
664
665
def effective_viscosity(self, dt: float) -> ufl.core.expr.Expr:
    """Effective viscosity used to impose boundary conditions on displacement
    weakly through a Nitsche penalty term in the viscosity_term of
    momentum_equation.py"""

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

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

Bases: QuasiCompressibleInternalVariableApproximation

Compressible viscoelastic rheology via the internal variable formulation.

N.b. this class includes the two additional terms neglected in the QuasiCompressibleInternalVariableApproximation: compressible buoyancy and the volume integral after integrating the advection of hydrostatic prestress term by parts.

For more information on the methodology and references please see 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 | 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
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 
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
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 
793
794
795
796
797
798
799
800
801
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)