Reproducing rheology of Morency and Doin, 2004

This section was contributed by Jonathan Perry-Houts

Modeling interactions between the upper mantle and the lithosphere can be difficult because of the dynamic range of temperatures and pressures involved. Many simple material models will assign very high viscosities at low temperature thermal boundary layers. The pseudo-brittle rheology described in Morency and Doin [2004] was developed to limit the strength of lithosphere at low temperature. The effective viscosity can be described as the harmonic mean of two non-Newtonian rheologies:

\[v_{\text{eff}} = \left(\frac{1}{v_{\text{eff}}^v}+\frac{1}{v_{\text{eff}}^p}\right)^{-1}\]

where

\[\begin{split}\begin{aligned} v_{\text{eff}}^v = B \left(\frac{\dot{\epsilon}}{\dot{\epsilon}_\text{ref}}\right)^{-1+1/n_v} \exp\left(\frac{E_a +V_a \rho_m g z}{n_v R T}\right), \\ v_{\text{eff}}^p = (\tau_0 + \gamma \rho_m g z) \left( \frac{\dot{\epsilon}^{-1+1/n_p}} {\dot{\epsilon}_\text{ref}^{1/n_p}} \right), \end{aligned}\end{split}\]

where \(B\) is a scaling constant; \(\dot{\epsilon}\) is defined as the quadratic sum of the second invariant of the strain rate tensor and a minimum strain rate, \(\dot{\epsilon}_0\); \(\dot{\epsilon}_\text{ref}\) is a reference strain rate; \(n_v\), and \(n_p\) are stress exponents; \(E_a\) is the activation energy; \(V_a\) is the activation volume; \(\rho_m\) is the mantle density; \(R\) is the gas constant; \(T\) is temperature; \(\tau_0\) is the cohesive strength of rocks at the surface; \(\gamma\) is a coefficient of yield stress increase with depth; and \(z\) is depth.

By limiting the strength of the lithosphere at low temperature, this rheology allows one to more realistically model processes like lithospheric delamination and foundering in the presence of weak crustal layers. A similar model setup to the one described in Morency and Doin [2004] can be reproduced with the files in the directory cookbooks/morency_doin_2004. In particular, the following sections of the input file are important to reproduce the setup:

Note

Morency and Doin [2004] defines the second invariant of the strain rate in a nonstandard way. The formulation in the paper is given as \(\epsilon_{II} = \sqrt{\frac{1}{2} (\epsilon_{11}^2 + \epsilon_{12}^2)}\), where \(\epsilon\) is the strain rate tensor. For consistency, that is also the formulation implemented in ASPECT. Because of this irregularity it is inadvisable to use this material model for purposes beyond reproducing published results.

Note

The viscosity profile in Figure 1 of Morency and Doin [2004] appears to be wrong. The published parameters do not reproduce these viscosities; it is unclear why. The values used here get very close. See Fig. 95 for an approximate reproduction of the original figure.

subsection Geometry model
  set Model name = box

  subsection Box
    set X extent      = 3000e3
    set Y extent      = 750e3
    set X repetitions = 4
  end
end

subsection Compositional fields
  set Number of fields = 2
  set Names of fields  = upper_crust, lower_crust
end

subsection Initial composition model
  set Model name = function

  subsection Function
    set Variable names = x,y
    set Function expression = if(y>=725e3,1,0);if((y<725e3&y>700e3),1,0)
  end
end

subsection Initial temperature model
  set Model name = function

  subsection Function
    set Variable names = x,y
    set Function constants = h=750e3, w=3000e3, mantleT=1350 # deg C
    set Function expression = \
      if( y < 100e3, \
        (100e3-y)/100e3*(1600-mantleT)+mantleT+293, \
        if(y>650e3, \
          (h-y)/(100e3)*mantleT+293, \
          mantleT+293))
  end
end

subsection Material model
  set Model name = Morency and Doin

  subsection Morency and Doin
    set Densities = 3300,2920,2920
    set Activation energies = 500,320,320
    set Coefficient of yield stress increase with depth = 0.25
    set Thermal expansivities = 3.5e-5
    set Stress exponents for viscous rheology = 3
    set Stress exponents for plastic rheology = 30
    set Thermal diffusivity = 0.8e-6
    set Heat capacity = 1.25e3
    set Activation volume = 6.4e-6
    set Reference strain rate = 6.4e-16
    set Preexponential constant for viscous rheology law = 7e11 ## Value used in paper is 1.24e14
    set Cohesive strength of rocks at the surface = 117
    set Reference temperature = 293
    set Minimum strain rate = 5e-19                             ## Value used in paper is 1.4e-20
  end
end

Fig. 95 Approximate reproduction of figure 1 from Morency and Doin [2004] using the ‘morency doin’ material model with almost all default parameters. Note the low-viscosity Moho, enabled by the low activation energy of the crustal component.