# A propagator matrix method for the Rayleigh–Taylor instability of multiple layers: a case study on crustal delamination in the early Earth

A propagator matrix method for the Rayleigh–Taylor instability of multiple layers: a case study... Summary The dispersion relation of the Rayleigh–Taylor instability, a gravitational instability associated with unstable density stratification, is of profound importance in various geophysical contexts. When more than two layers are involved, a semi-analytical technique based on the biharmonic formulation of Stokes flow has been extensively used to obtain such dispersion relation. However, this technique may become cumbersome when applied to lithospheric dynamics, where a number of layers are necessary to represent the continuous variation of viscosity over many orders of magnitude. Here, we present an alternative and more efficient method based on the propagator matrix formulation of Stokes flow. With this approach, the original instability problem is reduced to a compact eigenvalue equation whose size is solely determined by the number of primary density contrasts. We apply this new technique to the stability of the early crust, and combined with the Monte Carlo sensitivity analysis, we derive an empirical formula to compute the growth rate of the Rayleigh–Taylor instability for this particular geophysical setting. Our analysis indicates that the likelihood of crustal delamination hinges critically on the effective viscosity of eclogite. Numerical approximations and analysis, Dynamics of lithosphere and mantle, Rheology: crust and lithosphere 1 INTRODUCTION The Rayleigh–Taylor instability is a classical problem in fluid mechanics (Taylor 1950; Chandrasekhar 1961). This refers to the instability of an interface between two fluid layers with the upper layer being denser than the lower one. Given that this situation is always gravitationally unstable in the absence of surface tension, the problem amounts to finding how the growth rate of perturbations depends on their wavelengths or calculating the maximum growth rate and corresponding wavelength. The Rayleigh–Taylor instability is important in various geophysical contexts, such as salt dome formation (Selig 1965; Zaleski & Julien 1992), the delamination of crust or lithosphere (Fletcher & Hallet 1983; Bassi & Bonnin 1988; Houseman & Molnar 1997; Conrad & Molnar 1997), mantle plumes (Whitehead & Luther 1975) and the formation of the Earth's core (Ida et al. 1987; Stevenson 1990). The relation between the growth rate and wavelength of perturbations is well established for simple two-layer cases. For the delamination of crust or lithosphere, however, considering the instability of multiple layers becomes important, and it is common to use a semi-analytical technique developed by Bassi & Bonnin (1988). It is based on a set of fourth-order ordinary differential equations, thereby requiring to determine four unknown parameters for each layer. These unknown coefficients are then used to construct an eigenvalue problem, the solution of which contains the growth rate under consideration. When the number of layers is N, this method requires inverting a (4N + 2) × (4N + 2) matrix as well as solving an eigenvalue problem with a N × N matrix, the computational complexity of which are O(100N3) and O(N3), respectively. As a consequence, it becomes cumbersome to use this method to study the instability of lithosphere, where viscosity varies continuously over many orders of magnitude, because a number of layers are required to accurately handle such viscosity variations. In this paper, we present an alternative formulation using propagator matrices, which allows us to construct a compact eigenvalue problem; in this method the size of the matrix involved is determined by the number of primary density jumps, not the total number of layers. Propagator matrices have long been used in geodynamics to provide semi-analytical solutions for Stokes flow when viscosity varies only vertically (e.g. Hager & O'Connell 1981; Forte 2000). The computational efficiency of our approach is achieved by utilizing the fact that propagator matrices provide an analytically compact representation of fluid flow associated with the Rayleigh–Taylor instability of multiple layers. As a worked example, we apply this new approach to the problem of crustal delamination in the early Earth. The early oceanic crust is believed to be considerably thicker than the present oceanic crust as a consequence of high mantle potential temperature (∼1600 ° C, Herzberg et al. 2010). At the bottom part of such thick crust, basalt could transform to eclogite, which is denser than the underlying lithospheric mantle, thereby leading to the gravitational instability of the crust–mantle boundary (Johnson et al. 2014). The structure of the paper is following. We begin with the theoretical formulation of the Rayleigh–Taylor instability for a system with an arbitrary number of layers using the propagator matrix method for instantaneous Stokes flow. We then test our method with various analytical solutions. An application to the stability of the early crust is presented next. The Rayleigh–Taylor instability of multiple layers depends on a number of parameters, and it can be challenging to understand how the growth rate varies in the high-dimensional parameter space. In this worked example, therefore, we also discuss a practical approach to identify a set of the most important parameters by combining Monte Carlo sampling and linear regression. Our results suggest that the likelihood of crustal delamination in the early Earth depends critically on the viscosity of eclogite, which can change drastically along with the cooling of the early crust. We close with some important caveats for the problem of crustal delamination. 2 THEORETICAL FORMULATION The governing equations for an incompressible viscous fluid at the limit of infinite Prandtl number include the equation of continuity   $$\nabla \cdot {\bf V}=0,$$ (1)the equation of motion   $$0=\nabla \cdot \mathbf {\Sigma }+{\bf F},$$ (2)and the constitutive relation   $$\mathbf {\Sigma }=-P{\bf I}+\eta \left[\nabla \textbf {V}+(\nabla {\bf V})^{T}\right],$$ (3)where V is the total velocity field, Σ is the total stress, F is the body force, P is the total pressure, I is the identity matrix and η is the dynamic viscosity. To focus on the linear Rayleigh–Taylor instability, we express the flow-related entities as first-order perturbations from the hydrostatic reference:   \begin{eqnarray} {\bf V}&=&{\bf V}_{0}+\epsilon {\bf v}, \end{eqnarray} (4)  \begin{eqnarray} \mathbf {\Sigma }&=&\mathbf {\Sigma }_{0}+\epsilon \mathbf {\sigma }, \end{eqnarray} (5)  \begin{eqnarray} P&=&P_{0}+\epsilon p, \end{eqnarray} (6)  \begin{eqnarray} {\bf F}&=&{\bf F}_{0}+\epsilon {\bf f}, \end{eqnarray} (7)where [V0, Σ0, P0, F0]T and [v, σ, p, f]T represent the hydrostatic reference state (V0 is trivially zero) and the perturbed state, respectively. Subtracting the hydrostatic reference state leads to   \begin{eqnarray} v_{j,j}&=&0, \end{eqnarray} (8)  \begin{eqnarray} \sigma _{ij,j}+f_{i}&=&0, \end{eqnarray} (9)  \begin{eqnarray} \sigma _{ij}&=&-p\delta _{ij}+\eta \left(v_{i,j}+v_{j,i}\right), \end{eqnarray} (10)where fi arises solely due to density contrasts at the interfaces of two adjacent fluids. The system of eqs (8)–(10) describing the dynamics of the perturbed state is used to calculate the growth rate of initial density perturbations. In this study, our formulation is based on 2-D Cartesian geometry, but our results are applicable to 3-D Cartesian geometry as well; in the regime of linear instability, the growth rate depends on the magnitude of the wavenumber vector of perturbations, so the spatial pattern of perturbations is not important (Chandrasekhar 1961). In this 2-D Cartesian geometry, we expand field variables in the Fourier series assuming the horizontal invariance of material properties:   $$\xi (x,z)=\sum _{m=1}^{\infty }\left[\xi ^{1}_{m}(z)\cos (k_{m}x)+\xi ^{2}_{m}(z)\sin (k_{m}x)\right],$$ (11)where km is the horizontal wavenumber. Utilizing the orthogonality of the trigonometric functions, we may derive two sets of equations corresponding to the two Fourier coefficients (i.e. $$\xi ^{1}_{m}$$ and $$\xi ^{2}_{m}$$) from eqs (8)–(10). Since we are interested in the growth rate of the instability not the complete flow field, we may proceed with one of these two sets. Hereinafter, the positive z direction points vertically downward. Eqs (8)–(10), together with eq. (11), lead to   \begin{eqnarray} Dw&=&-ku, \end{eqnarray} (12)  \begin{eqnarray} Du&=&kw+\frac{1}{\eta }\sigma _{xz}, \end{eqnarray} (13)  \begin{eqnarray} D\sigma _{zz}&=&-k\sigma _{xz}-\Delta \rho g, \end{eqnarray} (14)  \begin{eqnarray} D\sigma _{xz}&=&k\sigma _{zz}+4\eta k^{2}u, \end{eqnarray} (15)where w, u, σzz and σxz are the Fourier coefficients of vertical velocity, horizontal velocity, vertical normal stress and shear stress, respectively, and Δρ represents density contrasts at fluid interfaces. The subscript m is omitted for simplicity. Eqs (12)–(15) describe the flow field induced by the density perturbation Δρ. Whereas these equations may be combined into a fourth-order ordinary differential equation (e.g. Bassi & Bonnin 1988), we retain the set of first-order differential equations to use the propagator matrix approach. First, the above Fourier coefficients are transformed to the following form for mathematical convenience:   \begin{eqnarray} y_{1}&=&w, \end{eqnarray} (16)  \begin{eqnarray} y_{2}&=&u, \end{eqnarray} (17)  \begin{eqnarray} y_{3}&=&\frac{\sigma _{zz}}{2\eta _{0}k}, \end{eqnarray} (18)  \begin{eqnarray} y_{4}&=&\frac{\sigma _{xz}}{2\eta _{0}k}, \end{eqnarray} (19)where η0 is the reference viscosity. Eqs (12)–(15) may then be written more compactly as   $$\frac{d\mathbf {y}}{dz}={\bf A}\mathbf {y}+\mathbf {b}.$$ (20)Here, y (≡ [y1, y2, y3, y4]T) is the velocity–stress vector, b (≡ [0, 0, Δρ g/2η0k, 0]T) is the buoyancy vector and   $$\mathbf {A}=\left[ \begin{array}{lccc} 0 &-k &\quad 0 & 0\\ k &0 &\quad 0 &\quad 2k/\eta ^{*}\\ 0 &0 &\quad0 &-k\\ 0 &\quad 2\eta ^{*}k &\quad k &0\\ \end{array} \right],$$ (21)where η* (≡ η/η0) is the non-dimensional viscosity. This equation has a simple exponential solution if A is a constant matrix. When viscosity is vertically varying, we may represent its variation by a stack of homogeneous layers, with each layer having a constant viscosity. Within each of such homogeneous layers, an exact solution to eq. (20) may be expressed as   $$\mathbf {y}(z)={\bf P}(z,z^{\prime })\mathbf {y}(z^{\prime })+\int _{z^{\prime }}^{z}{\bf P}(z,\xi ){\bf b}(\xi )d\xi ,$$ (22)where z΄ denotes the bottom of the layer, and P(z, z΄) is the propagator matrix expressed as   $${\bf P}(z,z^{\prime })=\sum _{n=0}^{\infty }\frac{(z-z^{\prime })^{n}}{n!}{\bf A}^{n}=\exp \left[{\bf A}(z-z^{\prime })\right].$$ (23) Let us now proceed to solve for the growth rate as a function of the wavenumber using the propagator matrices defined above (eqs 22 and 23). First, to construct the buoyancy vector b, we need to provide a mathematical description for density jumps arising at the interfaces. We assume that such density jumps may be represented by the Dirac delta functions as long as we confine ourselves to the regime of linear instability. In an M-layer system, we denote the amplitudes of interface deformations collectively as h (≡ [h1, h2, h3, ……, hM − 1]T). The growth or decay of these instabilities depends on the signs of the density jumps. The total density perturbation may be written as a sum over these individual jumps:   $$\Delta \rho =\sum _{i=1}^{M-1}h_{i}(\rho _{i}-\rho _{i+1})\delta (z-z_{i}).$$ (24) Now that we have defined the buoyancy vector, the next step entails the choice of the boundary condition. To make our derivation reasonably general, we assume that any two components of the field y are zero at a boundary; this can include both free-slip and rigid boundary conditions. Let us assume that lth and mth components are non-zero at the bottom boundary (i.e. $$y_M^l \ne 0$$ and $$y_M^m \ne 0$$), whereas ith and jth components are zero at the top boundary (i.e. $$y^{i}_{0}=0$$ and $$y^{j}_{0}=0$$). This boundary condition is utilized as follows. We propagate the velocity–stress vector y(zM) from the bottom boundary to the base of the (M − 1)th interface, that is, z = zM − 1 − ε, as:   $$\mathbf {y}(z_{M-1}-\epsilon )={\bf P}(z_{M-1}-\epsilon ,z_{M})\mathbf {y}(z_{M}).$$ (25)It may be observed that the discontinuity of density causes a jump in the normal stress across z = zM − 1. To account for this jump, we propagate the velocity–stress vector Y(zM − 1 − ε) across the interface, that is, from zM − 1 − ε to zM − 1 + ε:   $${\bf y}(z_{M-1}+\epsilon )={\bf P}(z_{M-1}+\epsilon ,z_{M-1}-\epsilon ){\bf y}(z_{M-1}-\epsilon )+\Bigg[0,0,\frac{h_{M-1}(\rho _{M-1}-\rho _{M})g}{2\eta ^{*}k},0\Bigg]^{T}.$$ (26)We may continue to propagate the velocity–stress vector in a similar manner until we reach the topmost boundary:   \begin{eqnarray} \mathbf {y}(z_{M-2}-\epsilon )&=&{\bf P}(z_{M-2}-\epsilon ,z_{M-1}+\epsilon )\mathbf {y}(z_{M-1}+\epsilon ),\nonumber \\ \mathbf {y}(z_{M-2}+\epsilon )&=&{\bf P}(z_{M-2}+\epsilon ,z_{M-2}-\epsilon )\mathbf {y}(z_{M-2}-\epsilon )+\Bigg[0,0,\frac{h_{M-2}(\rho _{M-2}-\rho _{M-1})g}{2\eta ^{*}k},0\Bigg]^{T},\nonumber \\ \vdots \nonumber \\ \mathbf {y}(z_{0})&=&{\bf P}(z_{0},z_{1}+\epsilon )\mathbf {y}(z_{1}+\epsilon ). \end{eqnarray} (27)The velocity–stress vector at the top boundary may be written in the following form using eqs (25)–(27) at the limit of ε → 0:   $$\mathbf {y}(z_{0})=\left[\prod _{i=1}^{M}{\bf P}(z_{i-1},z_{i})\right]\mathbf {y}(z_{M})+\sum _{n=1}^{M-1}\prod _{i=1}^{n}{\bf P}(z_{0},z_{i}){\bf b}_{i},$$ (28)where bi is the buoyancy vector arising at the ith interface, and we have used limε → 0P(zi − ε, zi + ε) = I. By applying the boundary condition to the velocity–stress vector expressed as eq. (28), we may derive two equations for the two unknown non-zero components of y(zM):   \begin{eqnarray} y^{l}_{M}P_{M}^{il}+y^{m}_{M}P_{M}^{im}&=&-\sum _{k=1}^{M-1}P_{k}^{i3}\frac{(\rho _{k}-\rho _{k+1})g}{2\eta _{0}k}h_{k}, \end{eqnarray} (29)  \begin{eqnarray} y^{l}_{M}P_{M}^{jl}+y^{m}_{M}P_{M}^{jm}&=&-\sum _{k=1}^{M-1}P_{k}^{j3}\frac{(\rho _{k}-\rho _{k+1})g}{2\eta _{0}k}h_{k}, \end{eqnarray} (30)where $$P^{ij}_{n}$$ refers to the (i, j) component of Pn, and $$\mathbf {P}_{n} \equiv \prod _{i=0}^{n-1}{\bf P}(z_{i},z_{i+1})={\bf P}(z_{0},z_{n})$$. Eqs (29) and (30) may be written in a more compact form as   $${\bf T}\mathbf {c}={\bf B}\mathbf {h},$$ (31)where $$\mathbf {c}=[y^{l}_{M},y^{m}_{M}]^{T}$$,   $${\bf T}= \left[\begin{array}{lc} P_{M}^{il} &\quad P_{M}^{im}\\ P_{M}^{jl} &\quad P_{M}^{jm} \end{array} \right],$$ (32)and   $${\bf B}= \left[\begin{array}{lccc} -P_{1}^{i3}\frac{(\rho _{1}-\rho _{2})g}{2\eta _{0}k} &\quad -P_{2}^{i3}\frac{(\rho _{2}-\rho _{3})g}{2\eta _{0}k} &\quad \dots &\quad -P_{M-1}^{i3}\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ -P_{1}^{j3}\frac{(\rho _{1}-\rho _{2})g}{2\eta _{0}k} &\quad -P_{2}^{j3}\frac{(\rho _{2}-\rho _{3})g}{2\eta _{0}k} &\quad \dots &\quad-P_{M-1}^{j3}\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ \end{array} \right].$$ (33)We may thus obtain c as   $$\mathbf {c}={\bf T}^{-1}{\bf B}\mathbf {h}.$$ (34) The vector c uniquely determines the velocity–stress vector at z = zM. Subsequently, we can obtain the velocity–stress vector at ∀z ∈ [z0, zM]. Using eqs (25)–(27), we may compute vertical velocity at every interface as   $$\mathbf {w}_{{\rm interface}}={\bf R}\mathbf {c}+{\bf B}_{I}\mathbf {h},$$ (35)where winterface = [w(z1), w(z2), …, w(zM − 1)]T,   $${\bf R}= \left[\begin{array}{lc} Q^{1l}_{1} &\quad Q^{1m}_{1}\\ Q^{1l}_{2} &\quad Q^{1m}_{2}\\ \vdots &\quad \vdots \\ Q^{1l}_{M-1} &\quad Q^{1m}_{M-1} \end{array} \right],$$ (36)  $${\bf B}_{I}= \left[\begin{array}{ccccc} 0 &\quad P^{13}(z_{1},z_{2})\frac{(\rho _{2}-\rho _{3})g}{2\eta _{0}k} &\quad P^{13}(z_{1},z_{3})\frac{(\rho _{3}-\rho _{4})g}{2\eta _{0}k} &\quad\dots &\quad P^{13}(z_{1},z_{M-1})\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ 0 &\quad 0 &\quad P^{13}(z_{2},z_{3})\frac{(\rho _{3}-\rho _{4})g}{2\eta _{0}k} &\quad\dots &\quad P^{13}(z_{2},z_{M-1})\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ \dots &\quad\dots &\quad\dots &\quad\dots &\quad\dots \\ 0 &\quad 0 &\quad0 &\quad\dots &\quad P^{13}(z_{M-2},z_{M-1})\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ 0 &\quad 0 &\quad0 &\quad\dots &\quad 0 \end{array} \right],$$ (37)and $${\bf Q}_{n} \equiv {\bf P}(z_{n},z_{M})=\prod _{i=n}^{M-1}{\bf P}(z_{i},z_{i+1})$$. The growth of the interface instability is equivalent to the temporal evolution of h, which may be related to the vertical velocity as   $${\bf w}_{{\rm interface}}=\frac{d \mathbf {h}}{dt}.$$ (38)By noting that h ∝ exp (qt) in the linear regime, eqs (34), (35) and (38) may be combined to form an eigenvalue problem as   $$({\bf R}{\bf T}^{-1}{\bf B}+{\bf B}_{I})\mathbf {h}=q \mathbf {h},$$ (39)where q is the growth rate. By solving this eigenvalue problem at a range of wavenumbers, we can construct the desired relation between the growth rate of instability and wavelength of perturbations. The size of the eigenvalue problem is defined by the number of primary density jumps, i.e. M − 1. The total number of layers used in our propagator matrix method may be much larger than M − 1 to represent large viscosity variations, but the size of the eigenvalue problem is determined only by the number of primary density jumps. The computational complexity of our semi-analytical technique is therefore O(M3), which is considerably less than O(N3), when N ≫ M. For geophysical applications, this efficiency may be an advantage, as viscosity varies continuously over many orders of magnitude. 3 SIMPLE EXAMPLES In this section, we provide some simple examples to show how our new method can handle different situations. First, we compare our semi-analytical solutions against exact solutions. Canright & Morris (1993) have derived the dispersion relation for the general two-layer case, when both layers have finite thicknesses and different viscosities. We have compared our solutions with such exact solutions. When one layer is infinitely thick, there also exist simpler solutions (e.g. Whitehead & Luther 1975; Conrad & Molnar 1997), and we also compare our solutions with such limiting cases. Then, we proceed to apply our method to a four-layer system that may represent a more realistic lithosphere model. With this four-layer example, we also consider the effect of phase transition. 3.1 Two-layer cases The velocity–stress vector at the top and bottom boundaries with the rigid boundary condition can be written as   \begin{eqnarray} \mathbf {y}(z_{0})=\bigg[0,0,y^{0}_{3},y^{0}_{4}\bigg]^{T}, \end{eqnarray} (40)  \begin{eqnarray} \mathbf {y}(z_{2})=\bigg[0,0,y^{2}_{3},y^{2}_{4}\bigg]^{T}. \end{eqnarray} (41)The thicknesses of the layers are d1 and d2, while their viscosities are η1 and η2, respectively. Before we proceed to solve the eigenvalue problem of eq. (39), it is convenient to non-dimensionalize the growth rate and the wavenumber as follows   \begin{eqnarray} q^{\prime }&=&q \left(\frac{(\rho _{1}-\rho _{2})gL}{2\eta _{0}}\right)^{-1}, \end{eqnarray} (42)  \begin{eqnarray} k^{\prime }&=&kL, \end{eqnarray} (43)where ρ1 and ρ2 are the densities of the upper and lower layers, respectively, η0 is the reference viscosity and it is set to be the viscosity of the top layer, and L is the characteristic length of the system, which is chosen to be min(d1, d2). First, we consider the case of d1 = d2 and compute the growth rate for a range of viscosity ratios (i.e. η2/η1). Similarly, we calculate the growth rates for a range of layer thickness ratios (i.e. d2/d1) when both layers have the same viscosity. These obtained growth rates are compared with the exact solutions (Figs 1a and b). The propagator matrix method can also handle the case of a thin layer overlying a semi-infinite half-space. Such a limiting case may be achieved approximately by letting d2 approach a reasonably high value such that d2 ≫ d1. Figs 1(c) and (e) show comparisons with the existing solution of Whitehead & Luther (1975), with d1 = 40 km and d2 = 800 km for two different viscosity combinations. Figure 1. View largeDownload slide Comparison of dispersion relations by propagator matrix (PM) method (shown by circles) against exact solutions (solid lines): (a) the case of two layers with the same thickness but different viscosities, (b) the case two layers with same viscosity but different thicknesses, (c) the case of a thin layer overlying a semi-infinite half-space (Whitehead & Luther 1975) for η2 ≥ η1, (d) the case of a slow linear decay of density within a thin layer overlying a semi-infinite half-space (Conrad & Molnar 1997) for η2 ≥ η1, (e) the case of a thin layer overlying a semi-infinite half-space (Whitehead & Luther 1975) for η2 ≤ η1 and (f) the case of a slow linear decay of density within a thin layer overlying a semi-infinite half-space (Conrad & Molnar 1997) for η2 ≤ η1. The chosen density profile is ρ(z) = ρ1 − (ρ1 − ρ2)z/D, where ρ1 = 400 kg m−3, ρ2 = 200 kg m−3 and D = 200 km is the thickness of the upper layer. In all cases, subscript 1 and 2 refer to the upper and lower layers, respectively. Figure 1. View largeDownload slide Comparison of dispersion relations by propagator matrix (PM) method (shown by circles) against exact solutions (solid lines): (a) the case of two layers with the same thickness but different viscosities, (b) the case two layers with same viscosity but different thicknesses, (c) the case of a thin layer overlying a semi-infinite half-space (Whitehead & Luther 1975) for η2 ≥ η1, (d) the case of a slow linear decay of density within a thin layer overlying a semi-infinite half-space (Conrad & Molnar 1997) for η2 ≥ η1, (e) the case of a thin layer overlying a semi-infinite half-space (Whitehead & Luther 1975) for η2 ≤ η1 and (f) the case of a slow linear decay of density within a thin layer overlying a semi-infinite half-space (Conrad & Molnar 1997) for η2 ≤ η1. The chosen density profile is ρ(z) = ρ1 − (ρ1 − ρ2)z/D, where ρ1 = 400 kg m−3, ρ2 = 200 kg m−3 and D = 200 km is the thickness of the upper layer. In all cases, subscript 1 and 2 refer to the upper and lower layers, respectively. In the examples considered so far, each layer has constant density and viscosity, but our method can be applied even when density varies continuously within each layer, and no primary jump occurs at the interface [e.g. Conrad & Molnar (1997)]. In such a case, we may approximate the density variation by a set of homogeneous layers with each having a constant density. The computed growth rate is compared with the analytical solution of Conrad & Molnar (1997) when density varies linearly within the upper layer (Figs 1d and f). 3.2 Four-layer cases We now move on to a four-layer case with each layer having finite thickness and different viscosity. Here, we assume the free-slip condition for both the top and bottom boundaries:   \begin{eqnarray} \mathbf {y}(z_{0})=\bigg[0,y^{0}_{2},y^{0}_{3},0\bigg]^{T}, \end{eqnarray} (44)  \begin{eqnarray} \mathbf {y}(z_{4})=\bigg[0,y^{4}_{2},y^{4}_{3},0\bigg]^{T}. \end{eqnarray} (45)The thicknesses, viscosities and densities of these four layers are denoted by di, ηi and ρi, respectively, where i ∈ [1, 4]. Similar to the two-layer case, we may choose the characteristic length scale to be mini(di). The density contrast at ith interface (i.e. ρi − ρi − 1) is denoted as Δρi. The reference viscosity (η0) is chosen to be meani(ηi). As an example, let us consider a case with (ρ1 − ρ2) < 0, (ρ2 − ρ3) > 0 and (ρ3 − ρ4) < 0. It is apparent from the signs of the density jumps that the first and the third interfaces (i.e. z = d1 and z = d1 + d2 + d3) stabilize themselves against any small perturbations to their flat reference states, whereas the interface at z = d1 + d2 is unstable to any small perturbations. The stable interfaces do not only resist self-deformation, they tend to suppress the instability arising at the unstable interface as a consequence of mass conservation. Now, the influence of a phase transition may be considered as follows. At the interface between the first and second layers (i.e. z = d1), the density jump is strictly negative (i.e. Δρ1 < 0). If the interface z = d1 corresponds to a phase boundary, and if the phase transition occurs fast enough compared to the interface deformation, any amount of a single phase crossing the phase boundary transforms to the other phase. The mathematical treatment of such physical situation is equivalent to Δρ1 = 0. This virtual zero density jump at z = d1 interface does not suppress the instability of z = d1 + d2 interface. In other words, the occurrence of phase transition at a stable interface enhances the growth of instability at the unstable interface, compared to the case of no phase transition. Comparison of these two situations is shown in the Fig. 2. Figure 2. View largeDownload slide (a) The growth rates computed in the presence (solid line) and absence (dashed line) of a phase boundary at the stable interface (z = d1) and (b) corresponding growth timescale. The parameters assumed are d1 = 35 km, d2 = 15 km, d3 = 100 km, d4 = 800 km, η1 = 1024 Pa s, η2 = 1023 Pa s, η3 = 1022 Pa s, η4 = 1020 Pa s, Δρ2 = 100 kg m−3 and Δρ3 = −10 kg m−3. Figure 2. View largeDownload slide (a) The growth rates computed in the presence (solid line) and absence (dashed line) of a phase boundary at the stable interface (z = d1) and (b) corresponding growth timescale. The parameters assumed are d1 = 35 km, d2 = 15 km, d3 = 100 km, d4 = 800 km, η1 = 1024 Pa s, η2 = 1023 Pa s, η3 = 1022 Pa s, η4 = 1020 Pa s, Δρ2 = 100 kg m−3 and Δρ3 = −10 kg m−3. 4 APPLICATION TO STABILITY OF THE ARCHEAN CRUST Oceanic crust in the Archean is likely to have been considerably thicker than the present-day oceanic crust owing to higher mantle potential temperature and thus higher degree of melting (e.g. Foley et al. 2003; Herzberg et al. 2010). The lowermost part of such thick crust may undergo metamorphic eclogitization, which increases its density and could subsequently lead to the gravitational instability of the crust–mantle boundary (Johnson et al. 2014). However, the actual physics of crustal delamination by such gravitational instability can be complex. The source of complexity lies mostly in uncertainties in mantle potential temperature, the thickness of the eclogite layer and the viscosity structure of the lithosphere. Because of these uncertainties, the model space of this multiple-layer Rayleigh–Taylor instability problem becomes large; we need to explore the consequences of each layer having different viscosities, thicknesses and densities. It can be challenging to understand the variation of the growth rate in a high-dimensional parameter space, but to access the geodynamical likelihood of crustal delamination, such consideration is essential. In this section, we provide a further application of our propagator matrix method to the problem of crustal delamination, and we also propose a practical approach to investigate efficiently the parameter dependence of the growth rate by combining Monte Carlo sampling and linear regression. Before presenting this new approach, we consider one particular example to set up a model for the Rayleigh–Taylor instability. 4.1 Thermal evolution and basalt–eclogite transition Almost every material property varies with temperature, so we first need to specify the thermal structure of our geodynamical model. The thermal evolution of oceanic lithosphere can be approximated by the standard half-space cooling model (Turcotte & Schubert 2002):   $$T(z,t)=T_{0}+(T_{m}-T_{0})\mbox{erf}\left(\frac{z}{2\sqrt{\kappa t}}\right)+z\left(\frac{\partial T}{\partial z}\right)_{S},$$ (46)where κ, T0 and Tm are the thermal diffusivity, the surface temperature and the initial mantle potential temperature, respectively. Our geodynamical model contains four layers: upper crust, lower crust, lithospheric mantle and asthenosphere. The lithospheric mantle and asthenosphere are considered 100 and 200 km thick, respectively. The whole crust is assumed to be 50 km thick, and the thickness of lower crust is determined by the basalt–eclogite transition, which is a function of pressure and temperature. Surface and mantle potential temperatures are assumed to be 0 ° C and 1650 ° C, respectively, along with a thermal diffusivity of 10−6 m2 s−1 and an adiabatic temperature gradient ((dT/dz)S) of 0.5 K km−1. Based on the eclogite stability field (Philpotts & Ague 2009), it can be observed that eclogite appears at the bottom part of the 50-km-thick crust only after ∼35 Myr (Fig. 3a). The thickness of the lower crust is time dependent, grows from 0 km at 35 Myr to ∼25 km at 100 Myr. Although this is admittedly a crude way to treat the basalt–eclogite transition (cf. Johnson et al. 2014), this degree of inaccuracy is unlikely to affect our main conclusion. Figure 3. View largeDownload slide (a) Lithospheric temperature profiles at 1, 5, 20, 50 and 100 Myr, assuming a half-space cooling model. The area with dashed boundary represents the approximate stability field for eclogite (Philpotts & Ague 2009). The dotted line represents the crust–mantle boundary. (b) Corresponding viscosity profiles. Eclogite appears only at 50 and 100 Myr. (c) The potential density profile. (d) The growth timescale is computed as a function of lithospheric age with three different viscosity profiles, i.e. continuous, average and minimum viscosity. Figure 3. View largeDownload slide (a) Lithospheric temperature profiles at 1, 5, 20, 50 and 100 Myr, assuming a half-space cooling model. The area with dashed boundary represents the approximate stability field for eclogite (Philpotts & Ague 2009). The dotted line represents the crust–mantle boundary. (b) Corresponding viscosity profiles. Eclogite appears only at 50 and 100 Myr. (c) The potential density profile. (d) The growth timescale is computed as a function of lithospheric age with three different viscosity profiles, i.e. continuous, average and minimum viscosity. We calculate temperature-dependent viscosity assuming dislocation creep, the flow law of which is expressed as   $$\dot{\varepsilon }=A\sigma ^{n}\exp \left(-\frac{E}{RT}\right),$$ (47)where $$\dot{\varepsilon }$$ is the strain rate, A is the pre-exponential factor, σ is the second invariant of the stress tensor, n is the stress exponent, E is the activation energy and R is the universal gas constant (e.g. Karato 2012, ch. 8). For the Rayleigh–Taylor instability, in which the stress grows from infinitesimal to some finite level, diffusion creep is most relevant. However, we lack experimental data of diffusion creep for garnet or eclogite, so we use the flow law of dislocation creep to estimate the effective viscosity, assuming a reference stress level of 0.01 MPa (Chu & Korenaga 2012). Also, we do not account for the pressure effect on viscosity, since only the shallow part of our model (<100 km) is important for crustal delamination. The rheology of mantle and upper crust is assumed to be represented by that of primary constituent minerals: olivine for mantle (Korenaga & Karato 2008) and plagioclase for upper crust (Rybacki et al. 2006). The rheology of eclogite is used for that of lower crust (Jin et al. 2001). Table 1 lists the flow-law parameters adopted in this study. Fig. 3(b) shows the viscosity profiles at different time instants. Table 1. The flow-law parameters assumed for upper crust (plagioclase), lower crust (eclogite) and mantle (olivine). Mineral or rock  E (kJ mol−1)  n  A (MPa−ns−1)  References  Plagioclase  641  3  1012.7  Rybacki et al. (2006)  Eclogite  480  3.4  103.3  Jin et al. (2001)  Olivine  610  4.94  106.09  Korenaga & Karato (2008)  Mineral or rock  E (kJ mol−1)  n  A (MPa−ns−1)  References  Plagioclase  641  3  1012.7  Rybacki et al. (2006)  Eclogite  480  3.4  103.3  Jin et al. (2001)  Olivine  610  4.94  106.09  Korenaga & Karato (2008)  View Large The evolution of density can be approximated by the linear relation   $$\Delta \rho (z,t)=\alpha \rho _{0}(z)[T(z,0)-T(z,t)],$$ (48)where α is thermal expansivity and ρ0 is the reference potential density (Korenaga & Korenaga 2016, appendix A) defined at initial potential temperature (Tm). The thermal expansivity is assumed to be 3 × 10−5 K−1. The reference potential density is set to 3000, 3340, 3290 and 3320 kg m−3, respectively, for upper crust, lower crust, lithospheric mantle and asthenosphere (Korenaga 2006; Korenaga & Korenaga 2016; Johnson et al. 2014). Note that we should use potential density instead of in situ density. When a material is isentropically brought down to a greater depth, for example, its temperature and density both increase, so when comparing densities at different depths, we always need to take into account this density variation along the adiabat, which can be achieved using potential density. We estimate the evolution of such potential density with the cooling of the lithosphere (Fig. 3c). As the lower crust becomes negatively buoyant with respect to the underlying lithospheric mantle after ∼35 Myr, only the cases of 50 and 100 Myr exhibit the unstable density stratification in Fig. 3(c). We use the average density for each layer obtained from the density profile while calculating the dispersion relation. This approach minimized the number of primary density contrasts; we also computed the dispersion relation with continuous density variation with no significant difference. The density jump at the upper–lower crust interface is set to zero to take into account the effects of phase transition (Section 3.2). Utilizing the propagator matrix technique, the growth timescale is computed using the assumed viscosity and density profiles for a range of lithospheric ages (Fig. 3d). It can be seen that the growth timescale is initially very long (∼8 Gyr) when the eclogitic lower crust first appears at ∼35 Myr; this is because the dense layer is very thin at the beginning. As the layer thickness increases with cooling, the growth timescale drops to ∼200 Myr at 40 Myr, and then gradually increases to 10 Gyr after 100 Myr of cooling. This increase originates in the temperature dependence of viscosity. In addition to the continuous viscosity profile, we also use the minimum and average viscosities of each layer, the result of which forms an envelope covering all the likely values of the growth timescale for the reasonable range of lithospheric viscosity. The minimum timescale of such envelope turns out to be ∼150 Myr corresponding to a lithospheric age of 40–50 Myr; the growth time is simply too long to cause crustal delamination, because the timescale increases to a few billion years at the age of 100 Myr. In other words, when the eclogitic crust becomes thick enough to be gravitationally unstable, it is already too stiff to flow by cooling, as a consequence of viscosity being strongly dependent on temperature. In this particular model, therefore, the delamination is unlikely to occur by the Rayleigh–Taylor instability. 4.2 Growth time formula for general cases In order to study the delamination of the early crust in a comprehensive manner, we must explore the plausible ranges of the model parameters involved and devise a practical method to compute the growth timescale. The following parameter ranges are deemed sufficiently wide for our purpose: d1 = 5–50 km, d2 = 5–50 km, d3 = 50–200 km, Δρ2 = 10–100 kg m−3 and Δρ3 = −30 to 0 kg m−3. The age of the oceanic lithosphere (t) is varied from 1 to 100 Myr, and the range of the reference stress (σref) is chosen between 0.001 and 1 MPa. With a chosen reference stress, the viscosity profile may be constructed following the flow law (eq. 47) at each time instant, which can be used to compute the minimum growth timescale and corresponding wavelength. First, we consider a reference case with d1 = 25 km, d2 = 25 km, d3 = 120 km, Δρ2 = 50 kg m−3, Δρ3 = −15 kg m−3, t = 50 Myr and log (σref) = −2. This particular set of input parameters results a minimum growth timescale of ∼300 Myr at a perturbation wavelength of ∼100 km. Now the model space to be explored can be defined around this reference state with the following standard deviation:   $$\sigma (\mathbf {x}_{1})=[10,10,35,15,10,35,1]^{T},$$ (49)where   $$\mathbf {x}_{1}=[\Delta z_{1},\Delta z_{2},\Delta z_{3},\Delta \rho _{2},\Delta \rho _{3},t,\log (\sigma _{{\rm ref}})]^{T}.$$ (50)Here, Δzi and Δρi are the thickness of the ith layer and the density jump at the ith interface, respectively. This model space is subjected to random sampling using the normal distribution, and a model vector (eq. 50) is constructed for each sample. This obtained set of random model vectors is then used to compute the dispersion relation. An ensemble of minimum growth timescales is obtained through 1000 such calculations. Even though we use the continuously varying viscosity while calculating the growth rate, we summarize a viscosity profile by calculating the log-mean viscosity of each layer, and we use such log-mean viscosities in the following linear regression. That is, instead of using the age of the oceanic crust and reference stress as two model parameters, we use a set of log-mean viscosities, so that we can directly measure the sensitivity of growth rate to viscosity. The log-mean viscosities corresponding to the reference state are as follows: η1 = 1024 Pa s, η2 = 1023 Pa s, η3 = 1020 Pa s and η4 = 1019 Pa s, respectively. The new model vector is defined with the log-mean viscosities ηi:   $$\mathbf {x}_{2}=[\Delta z_{1},\Delta z_{2},\Delta z_{3},\Delta \rho _{2},\Delta \rho _{3},\log (\eta _{1}),\log (\eta _{2}),\log (\eta _{3}),\log (\eta _{4})]^{T},$$ (51)with the standard deviation of   \begin{eqnarray} \sigma (\mathbf {x}_{2})&=&[10,10,35,15,10,2,1,0.5,0.5]^{T}. \end{eqnarray} (52)Despite the complexity of our model, the minimum growth timescale and corresponding wavelength of perturbations are found to be predicted with reasonable accuracy using linear functions:   \begin{eqnarray} \log _{10}(\tau )&=&a_{0}+\sum _{i=1}^{9}a_{i}X_{i}, \end{eqnarray} (53)  \begin{eqnarray} \log _{10}(\lambda )&=&b_{0}+\sum _{i=1}^{9}b_{i}X_{i}, \end{eqnarray} (54)where Xi is the mean subtracted input parameter that is normalized by the corresponding standard deviation (eq. 52). Constants a0 through b9 are estimated by the least-squares regression: a0 = 2.4350, b0 = 2.0407, a1 = 0.1197, a2 = 0.6143, a3 = 0.0093, a4 = 0.0897, a5 = 0.0058, a6 = 0.1099, a7 = 0.8303, a8 = 0.1319, a9 = 0.00042, and b1 = 0.0027, b2 = 0.0304, b3 = 0.0013, b4 = 0.0024, b5 = 0.0037, b6 = 0.0344, b7 = 0.086, b8 = 0.2115 and b9 = 0.00026. Fig. 4 demonstrates a reasonable correspondence between the predicted and actual values of the model outputs for the aforementioned 1000 simulations. This way of summarizing simulation results enables us not only to see the sensitivity of the growth timescale on the model parameters, but also to quickly reproduce the important results without doing the whole calculation. Keeping only the model parameters with high enough sensitivities, we may suggest the following expressions for the growth timescale and the wavelength of perturbations:   $$\log _{10}(\tau )= 2.435-0.061(\Delta z_{2}-25)-0.009(\Delta \rho _{2}-50)+0.830\left(\log _{10}(\eta _{2})-23\right) + 0.264\left(\log _{10}(\eta _{3})-20\right),$$ (55)  \begin{eqnarray} \log _{10}(\lambda )&=&2.041+0.003(\Delta z_{2}-25)-0.034\left(\log _{10}(\eta _{1})-24\right) \end{eqnarray} (56)  \begin{eqnarray} &&\quad \quad\quad \quad-0.086\left(\log _{10}(\eta _{2})-23\right)-0.423\left(\log _{10}(\eta _{3})-20\right), \end{eqnarray} (57)where Δz, Δρ, η, τ and λ are in km, kg m−3, Pa s, Myr and km, respectively. Figure 4. View largeDownload slide Correspondence between the predicted and actual values for (a) minimum growth timescale and (b) the wavelength of perturbations. Figure 4. View largeDownload slide Correspondence between the predicted and actual values for (a) minimum growth timescale and (b) the wavelength of perturbations. While the model space is high-dimensional, the growth timescale is primarily a function of the eclogite layer thickness, the density jump at crust–mantle interface, and most importantly, the viscosity of the eclogite layer. These formulae for growth timescale and wavelength of perturbation are helpful to assess the likelihood of crustal delamination in the early Earth. Based on thermomechanical calculations, Johnson et al. (2014) suggested the possibility of crustal delamination by the Rayleigh–Taylor instability with the timescale of a few million years. Although they have treated the thermodynamics of the basalt–eclogite transition in a detailed manner, they seemed to have neglected the impact of crustal rheology on the instability. As we have shown, the growth timescale varies strongly with the temperature-dependent viscosity of eclogite; cold and stiff lower crust becomes resistant to flow, and the growth timescale may be too high to allow the possibility of delamination. Note that eclogite is weaker than pure garnet aggregates (Karato et al. 1995; Jin et al. 2001), but even with the rheology of eclogite, the lower crust is too strong. It is possible, however, to reduce the effective viscosity of eclogite by invoking higher reference stress. The likelihood of crustal delamination thus seems to be limited to special tectonic settings. 5 DISCUSSION The propagator matrix formulation of the Rayleigh–Taylor instability has been tested and validated against analytical solutions. In this semi-analytical technique, the instability problem is reduced to a compact eigenvalue problem, whose size depends only on the number of primary density jumps. Therefore, the method can handle a range of different situations, including continuous variation of material properties. Also because of the problem being reduced to finding out the eigenvalue of a reasonably small matrix, the computational complexity of the problem (i.e. O(M3)) is considerably lower than that of the existing methods (e.g. Bassi & Bonnin 1988; Conrad & Molnar 1997). This computational efficiency enables us to execute a large number of computations quickly. Thus, we can explore the high-dimensional model space extensively to delineate the sensitivity of the growth rate on various model parameters. In other words, this matrix-based technique allows us to gain efficiently useful insight into the physics of the multiple-layer Rayleigh–Taylor instability in general. Though linear stability analysis does not provide any information on finite-amplitude phenomena, we can sweep the whole parameter space quickly and identify the model subspace that deserves further investigation. The propagator matrix formulation together with the Monte Carlo sensitivity analysis may be applied to other geophysical situations as well. For example, massive terrestrial planets with stagnant lid convection may feature a temperature profile that enters into the stability field of the eclogite after crust grows beyond a critical thickness (O'Rourke & Korenaga 2012). Thus, we may expect gravitational instability at the crust–mantle interface. Accessing the likelihood of crustal delamination in such cases through our analysis may lead to a strong conclusion whether these planets can escape the stagnant lid regime. Also, eclogite may be produced through pressure-induced phase transition in thick continental crust (Sobolev & Babeyko 2005), making it negatively buoyant with respect to the underlying mantle. The delamination of such thick continental crust (e.g. Bird 1979; Meissner & Mooney 1998; Jull & Kelemen 2001) may also be examined using our approach. Our analysis of early crust delamination underscores the importance of rheology to the dynamics of delamination. Therefore, it is of paramount importance to better understand lithospheric viscosity if we aim to be more conclusive about the likelihood of delamination by gravitational instabilities. With further deformation experiments combined with proper statistical analysis (e.g. Korenaga & Karato 2008; Mullet et al. 2015), the rheology of crust and mantle materials is expected to be constrained with smaller uncertainties. The Monte Carlo sensitivity analysis can point to the subset of material properties that are most important from the geodynamical perspective, thereby facilitating the feedback between mineral physics and geodynamics. ACKNOWLEDGEMENTS This study is based on work supported by the NASA through the NASA Astrobiology Institute under cooperative agreement no. NNA15BB03A issued through the Science Mission Directorate. We thank Gael Choblet and Neil Ribe for their careful reviews and constructive comments. REFERENCES Bassi G., Bonnin J., 1988. Rheological modelling and deformation instability of lithosphere under extension, Geophys. J. Int.  93 485– 504. Google Scholar CrossRef Search ADS   Bird P., 1979. Continental delamination and the Colorado Plateau, J. geophys. Res.  84 7561– 7571. Google Scholar CrossRef Search ADS   Canright D., Morris S., 1993. Buoyant instability of a viscous film over a passive fluid, J. Fluid Mech.  255 349– 372. Google Scholar CrossRef Search ADS   Chandrasekhar S., 1961. Hydrodynamic and Hydromagnetic Stability , Clarendon, Oxford. Chu X., Korenaga J., 2012. Olivine rheology, shear stress, and grain growth in the lithospheric mantle: geological constraints from the Kaapvaal craton, Earth. planet. Sci. Lett.  333 52– 62. Google Scholar CrossRef Search ADS   Conrad C.P., Molnar P., 1997. The growth of Rayleigh-Taylor type instabilities in the lithosphere for various rheological and density structures, Geophys. J. Int.  129 95– 112. Google Scholar CrossRef Search ADS   Fletcher R.C., Hallet B., 1983. Unstable extension of the lithosphere: a mechanical model for Basin and Range structure, J. geophys. Res.  88 7457– 7466. Google Scholar CrossRef Search ADS   Foley S.F., Buhre S., Jacob D.E., 2003. Evolution of the Archean crust by delamination and shallow subduction, Nature  421 249– 252. Google Scholar PubMed  Forte A.M., 2000. Seismic-geodynamic constraints on mantle flow: implications for layered convection, mantle viscosity, and seismic anisotropy in the deep mantle, in Earth's Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale , eds Karato S., Forte A., Liebermann R., Masters G., Stixrude L., Geophys. Monogr. Vol. 117, pp. 3– 36. AGU, Washington, DC. Google Scholar CrossRef Search ADS   Hager B.H., O'Connell R.J., 1981. A simple global model of plate dynamics and mantle convection, J. geophys. Res.  86 4843– 4867. Google Scholar CrossRef Search ADS   Herzberg C., Condie K., Korenaga J., 2010. Thermal history of the Earth and its petrological expression, Earth planet. Sci. Lett.  292 79– 88. Google Scholar CrossRef Search ADS   Houseman G.A., Molnar P., 1997. Gravitational (Rayleigh-Taylor) instability of a layer with non-linear viscosity and convective thinning of continental lithosphere, Geophys. J. Int.  128 125– 150. Google Scholar CrossRef Search ADS   Ida S., Nakagawa Y., Nakazawa K., 1987. The Earth's core formation due to the Rayleigh-Taylor instability, Icarus  69 239– 248. Google Scholar CrossRef Search ADS   Jin Z.-M., Zhang J., Green H., Jin S., 2001. Eclogite rheology: implications for subducted lithosphere, Geology  29 667– 670. Google Scholar CrossRef Search ADS   Johnson T.E., Brown M., Kaus B.J., VanTongeren J.A., 2014. Delamination and recycling of Archean crust caused by gravitational instabilities, Nature Geosci.  7 47– 52. Google Scholar CrossRef Search ADS   Jull M., Kelemen P.B., 2001. On the conditions for lower crustal convective instability, J. geophys. Res.  106 6423– 6446. Google Scholar CrossRef Search ADS   Karato S.-i., 2012. Deformation of Earth Materials: An Introduction to the Rheology of Solid Earth  Reprint edn, Cambridge University Press. Karato S.-i., Wang Z., Liu B., Fujino K., 1995. Plastic deformation of garnets: systematics and implications for the rheology of the mantle transition zone, Earth. planet. Sci. Lett.  130 13– 30. Google Scholar CrossRef Search ADS   Korenaga J., 2006. Archean geodynamics and the thermal evolution of Earth, in Archean Geodynamics and Environments  eds Benn K., Mareschal J-C., Condie K.C., Geophys. Monogr. Vol. 164, pp. 7– 32. AGU, Washington, DC. Google Scholar CrossRef Search ADS   Korenaga J., Karato S.-I., 2008. A new analysis of experimental data on olivine rheology, J. geophys. Res.  113 B02403, doi:10.1029/2007JB005100. Google Scholar CrossRef Search ADS   Korenaga T., Korenaga J., 2016. Evolution of young oceanic lithosphere and the meaning of seafloor subsidence rate, J. geophys. Res.  121 6315– 6332. Google Scholar CrossRef Search ADS   Meissner R., Mooney W., 1998. Weakness of the lower continental crust: a condition for delamination, uplift, and escape, Tectonophysics  296 47– 60. Google Scholar CrossRef Search ADS   Mullet B.G., Korenaga J., Karato S.-I., 2015. Markov chain Monte Carlo inversion for the rheology of olivine single crystals, J. geophys. Res.  120 3142– 3172. Google Scholar CrossRef Search ADS   O'Rourke J.G., Korenaga J., 2012. Terrestrial planet evolution in the stagnant-lid regime: size effects and the formation of self-destabilizing crust, Icarus  221 1043– 1060. Google Scholar CrossRef Search ADS   Philpotts A., Ague J., 2009. Principles of Igneous and Metamorphic Petrology  Cambridge University Press. Google Scholar CrossRef Search ADS   Rybacki E., Gottschalk M., Wirth R., Dresen G., 2006. Influence of water fugacity and activation volume on the flow properties of fine-grained anorthite aggregates, J. geophys. Res.  111 B03203, doi:10.1029/2005JB003663. Google Scholar CrossRef Search ADS   Selig F., 1965. A theoretical prediction of salt dome patterns, Geophysics  30 633– 643. Google Scholar CrossRef Search ADS   Sobolev S.V., Babeyko A.Y., 2005. What drives orogeny in the Andes?, Geology  33 617– 620. Google Scholar CrossRef Search ADS   Stevenson D., 1990. Fluid dynamics of core formation, in Origin of the Earth , pp. 231– 249, eds Newsom H., Jones J.E., Oxford University Press, London. Taylor G., 1950. The instability of liquid surfaces when accelerated in a direction perpendicular to their planes, Proc. R. Soc. Lond. A  201 192– 196. Google Scholar CrossRef Search ADS   Turcotte D., Schubert G., 2002. Geodynamics  Cambridge University Press. Google Scholar CrossRef Search ADS   Whitehead J.A., Luther D.S., 1975. Dynamics of laboratory diapir and plume models, J. geophys. Res.  80 705– 717. Google Scholar CrossRef Search ADS   Zaleski S., Julien P., 1992. Numerical simulation of Rayleigh-Taylor instability for single and multiple salt diapirs, Tectonophysics  206 55– 69. Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# A propagator matrix method for the Rayleigh–Taylor instability of multiple layers: a case study on crustal delamination in the early Earth

, Volume 212 (3) – Mar 1, 2018
12 pages

/lp/ou_press/a-propagator-matrix-method-for-the-rayleigh-taylor-instability-of-mZtoPr8GUd
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx513
Publisher site
See Article on Publisher Site

### Abstract

Summary The dispersion relation of the Rayleigh–Taylor instability, a gravitational instability associated with unstable density stratification, is of profound importance in various geophysical contexts. When more than two layers are involved, a semi-analytical technique based on the biharmonic formulation of Stokes flow has been extensively used to obtain such dispersion relation. However, this technique may become cumbersome when applied to lithospheric dynamics, where a number of layers are necessary to represent the continuous variation of viscosity over many orders of magnitude. Here, we present an alternative and more efficient method based on the propagator matrix formulation of Stokes flow. With this approach, the original instability problem is reduced to a compact eigenvalue equation whose size is solely determined by the number of primary density contrasts. We apply this new technique to the stability of the early crust, and combined with the Monte Carlo sensitivity analysis, we derive an empirical formula to compute the growth rate of the Rayleigh–Taylor instability for this particular geophysical setting. Our analysis indicates that the likelihood of crustal delamination hinges critically on the effective viscosity of eclogite. Numerical approximations and analysis, Dynamics of lithosphere and mantle, Rheology: crust and lithosphere 1 INTRODUCTION The Rayleigh–Taylor instability is a classical problem in fluid mechanics (Taylor 1950; Chandrasekhar 1961). This refers to the instability of an interface between two fluid layers with the upper layer being denser than the lower one. Given that this situation is always gravitationally unstable in the absence of surface tension, the problem amounts to finding how the growth rate of perturbations depends on their wavelengths or calculating the maximum growth rate and corresponding wavelength. The Rayleigh–Taylor instability is important in various geophysical contexts, such as salt dome formation (Selig 1965; Zaleski & Julien 1992), the delamination of crust or lithosphere (Fletcher & Hallet 1983; Bassi & Bonnin 1988; Houseman & Molnar 1997; Conrad & Molnar 1997), mantle plumes (Whitehead & Luther 1975) and the formation of the Earth's core (Ida et al. 1987; Stevenson 1990). The relation between the growth rate and wavelength of perturbations is well established for simple two-layer cases. For the delamination of crust or lithosphere, however, considering the instability of multiple layers becomes important, and it is common to use a semi-analytical technique developed by Bassi & Bonnin (1988). It is based on a set of fourth-order ordinary differential equations, thereby requiring to determine four unknown parameters for each layer. These unknown coefficients are then used to construct an eigenvalue problem, the solution of which contains the growth rate under consideration. When the number of layers is N, this method requires inverting a (4N + 2) × (4N + 2) matrix as well as solving an eigenvalue problem with a N × N matrix, the computational complexity of which are O(100N3) and O(N3), respectively. As a consequence, it becomes cumbersome to use this method to study the instability of lithosphere, where viscosity varies continuously over many orders of magnitude, because a number of layers are required to accurately handle such viscosity variations. In this paper, we present an alternative formulation using propagator matrices, which allows us to construct a compact eigenvalue problem; in this method the size of the matrix involved is determined by the number of primary density jumps, not the total number of layers. Propagator matrices have long been used in geodynamics to provide semi-analytical solutions for Stokes flow when viscosity varies only vertically (e.g. Hager & O'Connell 1981; Forte 2000). The computational efficiency of our approach is achieved by utilizing the fact that propagator matrices provide an analytically compact representation of fluid flow associated with the Rayleigh–Taylor instability of multiple layers. As a worked example, we apply this new approach to the problem of crustal delamination in the early Earth. The early oceanic crust is believed to be considerably thicker than the present oceanic crust as a consequence of high mantle potential temperature (∼1600 ° C, Herzberg et al. 2010). At the bottom part of such thick crust, basalt could transform to eclogite, which is denser than the underlying lithospheric mantle, thereby leading to the gravitational instability of the crust–mantle boundary (Johnson et al. 2014). The structure of the paper is following. We begin with the theoretical formulation of the Rayleigh–Taylor instability for a system with an arbitrary number of layers using the propagator matrix method for instantaneous Stokes flow. We then test our method with various analytical solutions. An application to the stability of the early crust is presented next. The Rayleigh–Taylor instability of multiple layers depends on a number of parameters, and it can be challenging to understand how the growth rate varies in the high-dimensional parameter space. In this worked example, therefore, we also discuss a practical approach to identify a set of the most important parameters by combining Monte Carlo sampling and linear regression. Our results suggest that the likelihood of crustal delamination in the early Earth depends critically on the viscosity of eclogite, which can change drastically along with the cooling of the early crust. We close with some important caveats for the problem of crustal delamination. 2 THEORETICAL FORMULATION The governing equations for an incompressible viscous fluid at the limit of infinite Prandtl number include the equation of continuity   $$\nabla \cdot {\bf V}=0,$$ (1)the equation of motion   $$0=\nabla \cdot \mathbf {\Sigma }+{\bf F},$$ (2)and the constitutive relation   $$\mathbf {\Sigma }=-P{\bf I}+\eta \left[\nabla \textbf {V}+(\nabla {\bf V})^{T}\right],$$ (3)where V is the total velocity field, Σ is the total stress, F is the body force, P is the total pressure, I is the identity matrix and η is the dynamic viscosity. To focus on the linear Rayleigh–Taylor instability, we express the flow-related entities as first-order perturbations from the hydrostatic reference:   \begin{eqnarray} {\bf V}&=&{\bf V}_{0}+\epsilon {\bf v}, \end{eqnarray} (4)  \begin{eqnarray} \mathbf {\Sigma }&=&\mathbf {\Sigma }_{0}+\epsilon \mathbf {\sigma }, \end{eqnarray} (5)  \begin{eqnarray} P&=&P_{0}+\epsilon p, \end{eqnarray} (6)  \begin{eqnarray} {\bf F}&=&{\bf F}_{0}+\epsilon {\bf f}, \end{eqnarray} (7)where [V0, Σ0, P0, F0]T and [v, σ, p, f]T represent the hydrostatic reference state (V0 is trivially zero) and the perturbed state, respectively. Subtracting the hydrostatic reference state leads to   \begin{eqnarray} v_{j,j}&=&0, \end{eqnarray} (8)  \begin{eqnarray} \sigma _{ij,j}+f_{i}&=&0, \end{eqnarray} (9)  \begin{eqnarray} \sigma _{ij}&=&-p\delta _{ij}+\eta \left(v_{i,j}+v_{j,i}\right), \end{eqnarray} (10)where fi arises solely due to density contrasts at the interfaces of two adjacent fluids. The system of eqs (8)–(10) describing the dynamics of the perturbed state is used to calculate the growth rate of initial density perturbations. In this study, our formulation is based on 2-D Cartesian geometry, but our results are applicable to 3-D Cartesian geometry as well; in the regime of linear instability, the growth rate depends on the magnitude of the wavenumber vector of perturbations, so the spatial pattern of perturbations is not important (Chandrasekhar 1961). In this 2-D Cartesian geometry, we expand field variables in the Fourier series assuming the horizontal invariance of material properties:   $$\xi (x,z)=\sum _{m=1}^{\infty }\left[\xi ^{1}_{m}(z)\cos (k_{m}x)+\xi ^{2}_{m}(z)\sin (k_{m}x)\right],$$ (11)where km is the horizontal wavenumber. Utilizing the orthogonality of the trigonometric functions, we may derive two sets of equations corresponding to the two Fourier coefficients (i.e. $$\xi ^{1}_{m}$$ and $$\xi ^{2}_{m}$$) from eqs (8)–(10). Since we are interested in the growth rate of the instability not the complete flow field, we may proceed with one of these two sets. Hereinafter, the positive z direction points vertically downward. Eqs (8)–(10), together with eq. (11), lead to   \begin{eqnarray} Dw&=&-ku, \end{eqnarray} (12)  \begin{eqnarray} Du&=&kw+\frac{1}{\eta }\sigma _{xz}, \end{eqnarray} (13)  \begin{eqnarray} D\sigma _{zz}&=&-k\sigma _{xz}-\Delta \rho g, \end{eqnarray} (14)  \begin{eqnarray} D\sigma _{xz}&=&k\sigma _{zz}+4\eta k^{2}u, \end{eqnarray} (15)where w, u, σzz and σxz are the Fourier coefficients of vertical velocity, horizontal velocity, vertical normal stress and shear stress, respectively, and Δρ represents density contrasts at fluid interfaces. The subscript m is omitted for simplicity. Eqs (12)–(15) describe the flow field induced by the density perturbation Δρ. Whereas these equations may be combined into a fourth-order ordinary differential equation (e.g. Bassi & Bonnin 1988), we retain the set of first-order differential equations to use the propagator matrix approach. First, the above Fourier coefficients are transformed to the following form for mathematical convenience:   \begin{eqnarray} y_{1}&=&w, \end{eqnarray} (16)  \begin{eqnarray} y_{2}&=&u, \end{eqnarray} (17)  \begin{eqnarray} y_{3}&=&\frac{\sigma _{zz}}{2\eta _{0}k}, \end{eqnarray} (18)  \begin{eqnarray} y_{4}&=&\frac{\sigma _{xz}}{2\eta _{0}k}, \end{eqnarray} (19)where η0 is the reference viscosity. Eqs (12)–(15) may then be written more compactly as   $$\frac{d\mathbf {y}}{dz}={\bf A}\mathbf {y}+\mathbf {b}.$$ (20)Here, y (≡ [y1, y2, y3, y4]T) is the velocity–stress vector, b (≡ [0, 0, Δρ g/2η0k, 0]T) is the buoyancy vector and   $$\mathbf {A}=\left[ \begin{array}{lccc} 0 &-k &\quad 0 & 0\\ k &0 &\quad 0 &\quad 2k/\eta ^{*}\\ 0 &0 &\quad0 &-k\\ 0 &\quad 2\eta ^{*}k &\quad k &0\\ \end{array} \right],$$ (21)where η* (≡ η/η0) is the non-dimensional viscosity. This equation has a simple exponential solution if A is a constant matrix. When viscosity is vertically varying, we may represent its variation by a stack of homogeneous layers, with each layer having a constant viscosity. Within each of such homogeneous layers, an exact solution to eq. (20) may be expressed as   $$\mathbf {y}(z)={\bf P}(z,z^{\prime })\mathbf {y}(z^{\prime })+\int _{z^{\prime }}^{z}{\bf P}(z,\xi ){\bf b}(\xi )d\xi ,$$ (22)where z΄ denotes the bottom of the layer, and P(z, z΄) is the propagator matrix expressed as   $${\bf P}(z,z^{\prime })=\sum _{n=0}^{\infty }\frac{(z-z^{\prime })^{n}}{n!}{\bf A}^{n}=\exp \left[{\bf A}(z-z^{\prime })\right].$$ (23) Let us now proceed to solve for the growth rate as a function of the wavenumber using the propagator matrices defined above (eqs 22 and 23). First, to construct the buoyancy vector b, we need to provide a mathematical description for density jumps arising at the interfaces. We assume that such density jumps may be represented by the Dirac delta functions as long as we confine ourselves to the regime of linear instability. In an M-layer system, we denote the amplitudes of interface deformations collectively as h (≡ [h1, h2, h3, ……, hM − 1]T). The growth or decay of these instabilities depends on the signs of the density jumps. The total density perturbation may be written as a sum over these individual jumps:   $$\Delta \rho =\sum _{i=1}^{M-1}h_{i}(\rho _{i}-\rho _{i+1})\delta (z-z_{i}).$$ (24) Now that we have defined the buoyancy vector, the next step entails the choice of the boundary condition. To make our derivation reasonably general, we assume that any two components of the field y are zero at a boundary; this can include both free-slip and rigid boundary conditions. Let us assume that lth and mth components are non-zero at the bottom boundary (i.e. $$y_M^l \ne 0$$ and $$y_M^m \ne 0$$), whereas ith and jth components are zero at the top boundary (i.e. $$y^{i}_{0}=0$$ and $$y^{j}_{0}=0$$). This boundary condition is utilized as follows. We propagate the velocity–stress vector y(zM) from the bottom boundary to the base of the (M − 1)th interface, that is, z = zM − 1 − ε, as:   $$\mathbf {y}(z_{M-1}-\epsilon )={\bf P}(z_{M-1}-\epsilon ,z_{M})\mathbf {y}(z_{M}).$$ (25)It may be observed that the discontinuity of density causes a jump in the normal stress across z = zM − 1. To account for this jump, we propagate the velocity–stress vector Y(zM − 1 − ε) across the interface, that is, from zM − 1 − ε to zM − 1 + ε:   $${\bf y}(z_{M-1}+\epsilon )={\bf P}(z_{M-1}+\epsilon ,z_{M-1}-\epsilon ){\bf y}(z_{M-1}-\epsilon )+\Bigg[0,0,\frac{h_{M-1}(\rho _{M-1}-\rho _{M})g}{2\eta ^{*}k},0\Bigg]^{T}.$$ (26)We may continue to propagate the velocity–stress vector in a similar manner until we reach the topmost boundary:   \begin{eqnarray} \mathbf {y}(z_{M-2}-\epsilon )&=&{\bf P}(z_{M-2}-\epsilon ,z_{M-1}+\epsilon )\mathbf {y}(z_{M-1}+\epsilon ),\nonumber \\ \mathbf {y}(z_{M-2}+\epsilon )&=&{\bf P}(z_{M-2}+\epsilon ,z_{M-2}-\epsilon )\mathbf {y}(z_{M-2}-\epsilon )+\Bigg[0,0,\frac{h_{M-2}(\rho _{M-2}-\rho _{M-1})g}{2\eta ^{*}k},0\Bigg]^{T},\nonumber \\ \vdots \nonumber \\ \mathbf {y}(z_{0})&=&{\bf P}(z_{0},z_{1}+\epsilon )\mathbf {y}(z_{1}+\epsilon ). \end{eqnarray} (27)The velocity–stress vector at the top boundary may be written in the following form using eqs (25)–(27) at the limit of ε → 0:   $$\mathbf {y}(z_{0})=\left[\prod _{i=1}^{M}{\bf P}(z_{i-1},z_{i})\right]\mathbf {y}(z_{M})+\sum _{n=1}^{M-1}\prod _{i=1}^{n}{\bf P}(z_{0},z_{i}){\bf b}_{i},$$ (28)where bi is the buoyancy vector arising at the ith interface, and we have used limε → 0P(zi − ε, zi + ε) = I. By applying the boundary condition to the velocity–stress vector expressed as eq. (28), we may derive two equations for the two unknown non-zero components of y(zM):   \begin{eqnarray} y^{l}_{M}P_{M}^{il}+y^{m}_{M}P_{M}^{im}&=&-\sum _{k=1}^{M-1}P_{k}^{i3}\frac{(\rho _{k}-\rho _{k+1})g}{2\eta _{0}k}h_{k}, \end{eqnarray} (29)  \begin{eqnarray} y^{l}_{M}P_{M}^{jl}+y^{m}_{M}P_{M}^{jm}&=&-\sum _{k=1}^{M-1}P_{k}^{j3}\frac{(\rho _{k}-\rho _{k+1})g}{2\eta _{0}k}h_{k}, \end{eqnarray} (30)where $$P^{ij}_{n}$$ refers to the (i, j) component of Pn, and $$\mathbf {P}_{n} \equiv \prod _{i=0}^{n-1}{\bf P}(z_{i},z_{i+1})={\bf P}(z_{0},z_{n})$$. Eqs (29) and (30) may be written in a more compact form as   $${\bf T}\mathbf {c}={\bf B}\mathbf {h},$$ (31)where $$\mathbf {c}=[y^{l}_{M},y^{m}_{M}]^{T}$$,   $${\bf T}= \left[\begin{array}{lc} P_{M}^{il} &\quad P_{M}^{im}\\ P_{M}^{jl} &\quad P_{M}^{jm} \end{array} \right],$$ (32)and   $${\bf B}= \left[\begin{array}{lccc} -P_{1}^{i3}\frac{(\rho _{1}-\rho _{2})g}{2\eta _{0}k} &\quad -P_{2}^{i3}\frac{(\rho _{2}-\rho _{3})g}{2\eta _{0}k} &\quad \dots &\quad -P_{M-1}^{i3}\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ -P_{1}^{j3}\frac{(\rho _{1}-\rho _{2})g}{2\eta _{0}k} &\quad -P_{2}^{j3}\frac{(\rho _{2}-\rho _{3})g}{2\eta _{0}k} &\quad \dots &\quad-P_{M-1}^{j3}\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ \end{array} \right].$$ (33)We may thus obtain c as   $$\mathbf {c}={\bf T}^{-1}{\bf B}\mathbf {h}.$$ (34) The vector c uniquely determines the velocity–stress vector at z = zM. Subsequently, we can obtain the velocity–stress vector at ∀z ∈ [z0, zM]. Using eqs (25)–(27), we may compute vertical velocity at every interface as   $$\mathbf {w}_{{\rm interface}}={\bf R}\mathbf {c}+{\bf B}_{I}\mathbf {h},$$ (35)where winterface = [w(z1), w(z2), …, w(zM − 1)]T,   $${\bf R}= \left[\begin{array}{lc} Q^{1l}_{1} &\quad Q^{1m}_{1}\\ Q^{1l}_{2} &\quad Q^{1m}_{2}\\ \vdots &\quad \vdots \\ Q^{1l}_{M-1} &\quad Q^{1m}_{M-1} \end{array} \right],$$ (36)  $${\bf B}_{I}= \left[\begin{array}{ccccc} 0 &\quad P^{13}(z_{1},z_{2})\frac{(\rho _{2}-\rho _{3})g}{2\eta _{0}k} &\quad P^{13}(z_{1},z_{3})\frac{(\rho _{3}-\rho _{4})g}{2\eta _{0}k} &\quad\dots &\quad P^{13}(z_{1},z_{M-1})\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ 0 &\quad 0 &\quad P^{13}(z_{2},z_{3})\frac{(\rho _{3}-\rho _{4})g}{2\eta _{0}k} &\quad\dots &\quad P^{13}(z_{2},z_{M-1})\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ \dots &\quad\dots &\quad\dots &\quad\dots &\quad\dots \\ 0 &\quad 0 &\quad0 &\quad\dots &\quad P^{13}(z_{M-2},z_{M-1})\frac{(\rho _{M-1}-\rho _{M})g}{2\eta _{0}k}\\ 0 &\quad 0 &\quad0 &\quad\dots &\quad 0 \end{array} \right],$$ (37)and $${\bf Q}_{n} \equiv {\bf P}(z_{n},z_{M})=\prod _{i=n}^{M-1}{\bf P}(z_{i},z_{i+1})$$. The growth of the interface instability is equivalent to the temporal evolution of h, which may be related to the vertical velocity as   $${\bf w}_{{\rm interface}}=\frac{d \mathbf {h}}{dt}.$$ (38)By noting that h ∝ exp (qt) in the linear regime, eqs (34), (35) and (38) may be combined to form an eigenvalue problem as   $$({\bf R}{\bf T}^{-1}{\bf B}+{\bf B}_{I})\mathbf {h}=q \mathbf {h},$$ (39)where q is the growth rate. By solving this eigenvalue problem at a range of wavenumbers, we can construct the desired relation between the growth rate of instability and wavelength of perturbations. The size of the eigenvalue problem is defined by the number of primary density jumps, i.e. M − 1. The total number of layers used in our propagator matrix method may be much larger than M − 1 to represent large viscosity variations, but the size of the eigenvalue problem is determined only by the number of primary density jumps. The computational complexity of our semi-analytical technique is therefore O(M3), which is considerably less than O(N3), when N ≫ M. For geophysical applications, this efficiency may be an advantage, as viscosity varies continuously over many orders of magnitude. 3 SIMPLE EXAMPLES In this section, we provide some simple examples to show how our new method can handle different situations. First, we compare our semi-analytical solutions against exact solutions. Canright & Morris (1993) have derived the dispersion relation for the general two-layer case, when both layers have finite thicknesses and different viscosities. We have compared our solutions with such exact solutions. When one layer is infinitely thick, there also exist simpler solutions (e.g. Whitehead & Luther 1975; Conrad & Molnar 1997), and we also compare our solutions with such limiting cases. Then, we proceed to apply our method to a four-layer system that may represent a more realistic lithosphere model. With this four-layer example, we also consider the effect of phase transition. 3.1 Two-layer cases The velocity–stress vector at the top and bottom boundaries with the rigid boundary condition can be written as   \begin{eqnarray} \mathbf {y}(z_{0})=\bigg[0,0,y^{0}_{3},y^{0}_{4}\bigg]^{T}, \end{eqnarray} (40)  \begin{eqnarray} \mathbf {y}(z_{2})=\bigg[0,0,y^{2}_{3},y^{2}_{4}\bigg]^{T}. \end{eqnarray} (41)The thicknesses of the layers are d1 and d2, while their viscosities are η1 and η2, respectively. Before we proceed to solve the eigenvalue problem of eq. (39), it is convenient to non-dimensionalize the growth rate and the wavenumber as follows   \begin{eqnarray} q^{\prime }&=&q \left(\frac{(\rho _{1}-\rho _{2})gL}{2\eta _{0}}\right)^{-1}, \end{eqnarray} (42)  \begin{eqnarray} k^{\prime }&=&kL, \end{eqnarray} (43)where ρ1 and ρ2 are the densities of the upper and lower layers, respectively, η0 is the reference viscosity and it is set to be the viscosity of the top layer, and L is the characteristic length of the system, which is chosen to be min(d1, d2). First, we consider the case of d1 = d2 and compute the growth rate for a range of viscosity ratios (i.e. η2/η1). Similarly, we calculate the growth rates for a range of layer thickness ratios (i.e. d2/d1) when both layers have the same viscosity. These obtained growth rates are compared with the exact solutions (Figs 1a and b). The propagator matrix method can also handle the case of a thin layer overlying a semi-infinite half-space. Such a limiting case may be achieved approximately by letting d2 approach a reasonably high value such that d2 ≫ d1. Figs 1(c) and (e) show comparisons with the existing solution of Whitehead & Luther (1975), with d1 = 40 km and d2 = 800 km for two different viscosity combinations. Figure 1. View largeDownload slide Comparison of dispersion relations by propagator matrix (PM) method (shown by circles) against exact solutions (solid lines): (a) the case of two layers with the same thickness but different viscosities, (b) the case two layers with same viscosity but different thicknesses, (c) the case of a thin layer overlying a semi-infinite half-space (Whitehead & Luther 1975) for η2 ≥ η1, (d) the case of a slow linear decay of density within a thin layer overlying a semi-infinite half-space (Conrad & Molnar 1997) for η2 ≥ η1, (e) the case of a thin layer overlying a semi-infinite half-space (Whitehead & Luther 1975) for η2 ≤ η1 and (f) the case of a slow linear decay of density within a thin layer overlying a semi-infinite half-space (Conrad & Molnar 1997) for η2 ≤ η1. The chosen density profile is ρ(z) = ρ1 − (ρ1 − ρ2)z/D, where ρ1 = 400 kg m−3, ρ2 = 200 kg m−3 and D = 200 km is the thickness of the upper layer. In all cases, subscript 1 and 2 refer to the upper and lower layers, respectively. Figure 1. View largeDownload slide Comparison of dispersion relations by propagator matrix (PM) method (shown by circles) against exact solutions (solid lines): (a) the case of two layers with the same thickness but different viscosities, (b) the case two layers with same viscosity but different thicknesses, (c) the case of a thin layer overlying a semi-infinite half-space (Whitehead & Luther 1975) for η2 ≥ η1, (d) the case of a slow linear decay of density within a thin layer overlying a semi-infinite half-space (Conrad & Molnar 1997) for η2 ≥ η1, (e) the case of a thin layer overlying a semi-infinite half-space (Whitehead & Luther 1975) for η2 ≤ η1 and (f) the case of a slow linear decay of density within a thin layer overlying a semi-infinite half-space (Conrad & Molnar 1997) for η2 ≤ η1. The chosen density profile is ρ(z) = ρ1 − (ρ1 − ρ2)z/D, where ρ1 = 400 kg m−3, ρ2 = 200 kg m−3 and D = 200 km is the thickness of the upper layer. In all cases, subscript 1 and 2 refer to the upper and lower layers, respectively. In the examples considered so far, each layer has constant density and viscosity, but our method can be applied even when density varies continuously within each layer, and no primary jump occurs at the interface [e.g. Conrad & Molnar (1997)]. In such a case, we may approximate the density variation by a set of homogeneous layers with each having a constant density. The computed growth rate is compared with the analytical solution of Conrad & Molnar (1997) when density varies linearly within the upper layer (Figs 1d and f). 3.2 Four-layer cases We now move on to a four-layer case with each layer having finite thickness and different viscosity. Here, we assume the free-slip condition for both the top and bottom boundaries:   \begin{eqnarray} \mathbf {y}(z_{0})=\bigg[0,y^{0}_{2},y^{0}_{3},0\bigg]^{T}, \end{eqnarray} (44)  \begin{eqnarray} \mathbf {y}(z_{4})=\bigg[0,y^{4}_{2},y^{4}_{3},0\bigg]^{T}. \end{eqnarray} (45)The thicknesses, viscosities and densities of these four layers are denoted by di, ηi and ρi, respectively, where i ∈ [1, 4]. Similar to the two-layer case, we may choose the characteristic length scale to be mini(di). The density contrast at ith interface (i.e. ρi − ρi − 1) is denoted as Δρi. The reference viscosity (η0) is chosen to be meani(ηi). As an example, let us consider a case with (ρ1 − ρ2) < 0, (ρ2 − ρ3) > 0 and (ρ3 − ρ4) < 0. It is apparent from the signs of the density jumps that the first and the third interfaces (i.e. z = d1 and z = d1 + d2 + d3) stabilize themselves against any small perturbations to their flat reference states, whereas the interface at z = d1 + d2 is unstable to any small perturbations. The stable interfaces do not only resist self-deformation, they tend to suppress the instability arising at the unstable interface as a consequence of mass conservation. Now, the influence of a phase transition may be considered as follows. At the interface between the first and second layers (i.e. z = d1), the density jump is strictly negative (i.e. Δρ1 < 0). If the interface z = d1 corresponds to a phase boundary, and if the phase transition occurs fast enough compared to the interface deformation, any amount of a single phase crossing the phase boundary transforms to the other phase. The mathematical treatment of such physical situation is equivalent to Δρ1 = 0. This virtual zero density jump at z = d1 interface does not suppress the instability of z = d1 + d2 interface. In other words, the occurrence of phase transition at a stable interface enhances the growth of instability at the unstable interface, compared to the case of no phase transition. Comparison of these two situations is shown in the Fig. 2. Figure 2. View largeDownload slide (a) The growth rates computed in the presence (solid line) and absence (dashed line) of a phase boundary at the stable interface (z = d1) and (b) corresponding growth timescale. The parameters assumed are d1 = 35 km, d2 = 15 km, d3 = 100 km, d4 = 800 km, η1 = 1024 Pa s, η2 = 1023 Pa s, η3 = 1022 Pa s, η4 = 1020 Pa s, Δρ2 = 100 kg m−3 and Δρ3 = −10 kg m−3. Figure 2. View largeDownload slide (a) The growth rates computed in the presence (solid line) and absence (dashed line) of a phase boundary at the stable interface (z = d1) and (b) corresponding growth timescale. The parameters assumed are d1 = 35 km, d2 = 15 km, d3 = 100 km, d4 = 800 km, η1 = 1024 Pa s, η2 = 1023 Pa s, η3 = 1022 Pa s, η4 = 1020 Pa s, Δρ2 = 100 kg m−3 and Δρ3 = −10 kg m−3. 4 APPLICATION TO STABILITY OF THE ARCHEAN CRUST Oceanic crust in the Archean is likely to have been considerably thicker than the present-day oceanic crust owing to higher mantle potential temperature and thus higher degree of melting (e.g. Foley et al. 2003; Herzberg et al. 2010). The lowermost part of such thick crust may undergo metamorphic eclogitization, which increases its density and could subsequently lead to the gravitational instability of the crust–mantle boundary (Johnson et al. 2014). However, the actual physics of crustal delamination by such gravitational instability can be complex. The source of complexity lies mostly in uncertainties in mantle potential temperature, the thickness of the eclogite layer and the viscosity structure of the lithosphere. Because of these uncertainties, the model space of this multiple-layer Rayleigh–Taylor instability problem becomes large; we need to explore the consequences of each layer having different viscosities, thicknesses and densities. It can be challenging to understand the variation of the growth rate in a high-dimensional parameter space, but to access the geodynamical likelihood of crustal delamination, such consideration is essential. In this section, we provide a further application of our propagator matrix method to the problem of crustal delamination, and we also propose a practical approach to investigate efficiently the parameter dependence of the growth rate by combining Monte Carlo sampling and linear regression. Before presenting this new approach, we consider one particular example to set up a model for the Rayleigh–Taylor instability. 4.1 Thermal evolution and basalt–eclogite transition Almost every material property varies with temperature, so we first need to specify the thermal structure of our geodynamical model. The thermal evolution of oceanic lithosphere can be approximated by the standard half-space cooling model (Turcotte & Schubert 2002):   $$T(z,t)=T_{0}+(T_{m}-T_{0})\mbox{erf}\left(\frac{z}{2\sqrt{\kappa t}}\right)+z\left(\frac{\partial T}{\partial z}\right)_{S},$$ (46)where κ, T0 and Tm are the thermal diffusivity, the surface temperature and the initial mantle potential temperature, respectively. Our geodynamical model contains four layers: upper crust, lower crust, lithospheric mantle and asthenosphere. The lithospheric mantle and asthenosphere are considered 100 and 200 km thick, respectively. The whole crust is assumed to be 50 km thick, and the thickness of lower crust is determined by the basalt–eclogite transition, which is a function of pressure and temperature. Surface and mantle potential temperatures are assumed to be 0 ° C and 1650 ° C, respectively, along with a thermal diffusivity of 10−6 m2 s−1 and an adiabatic temperature gradient ((dT/dz)S) of 0.5 K km−1. Based on the eclogite stability field (Philpotts & Ague 2009), it can be observed that eclogite appears at the bottom part of the 50-km-thick crust only after ∼35 Myr (Fig. 3a). The thickness of the lower crust is time dependent, grows from 0 km at 35 Myr to ∼25 km at 100 Myr. Although this is admittedly a crude way to treat the basalt–eclogite transition (cf. Johnson et al. 2014), this degree of inaccuracy is unlikely to affect our main conclusion. Figure 3. View largeDownload slide (a) Lithospheric temperature profiles at 1, 5, 20, 50 and 100 Myr, assuming a half-space cooling model. The area with dashed boundary represents the approximate stability field for eclogite (Philpotts & Ague 2009). The dotted line represents the crust–mantle boundary. (b) Corresponding viscosity profiles. Eclogite appears only at 50 and 100 Myr. (c) The potential density profile. (d) The growth timescale is computed as a function of lithospheric age with three different viscosity profiles, i.e. continuous, average and minimum viscosity. Figure 3. View largeDownload slide (a) Lithospheric temperature profiles at 1, 5, 20, 50 and 100 Myr, assuming a half-space cooling model. The area with dashed boundary represents the approximate stability field for eclogite (Philpotts & Ague 2009). The dotted line represents the crust–mantle boundary. (b) Corresponding viscosity profiles. Eclogite appears only at 50 and 100 Myr. (c) The potential density profile. (d) The growth timescale is computed as a function of lithospheric age with three different viscosity profiles, i.e. continuous, average and minimum viscosity. We calculate temperature-dependent viscosity assuming dislocation creep, the flow law of which is expressed as   $$\dot{\varepsilon }=A\sigma ^{n}\exp \left(-\frac{E}{RT}\right),$$ (47)where $$\dot{\varepsilon }$$ is the strain rate, A is the pre-exponential factor, σ is the second invariant of the stress tensor, n is the stress exponent, E is the activation energy and R is the universal gas constant (e.g. Karato 2012, ch. 8). For the Rayleigh–Taylor instability, in which the stress grows from infinitesimal to some finite level, diffusion creep is most relevant. However, we lack experimental data of diffusion creep for garnet or eclogite, so we use the flow law of dislocation creep to estimate the effective viscosity, assuming a reference stress level of 0.01 MPa (Chu & Korenaga 2012). Also, we do not account for the pressure effect on viscosity, since only the shallow part of our model (<100 km) is important for crustal delamination. The rheology of mantle and upper crust is assumed to be represented by that of primary constituent minerals: olivine for mantle (Korenaga & Karato 2008) and plagioclase for upper crust (Rybacki et al. 2006). The rheology of eclogite is used for that of lower crust (Jin et al. 2001). Table 1 lists the flow-law parameters adopted in this study. Fig. 3(b) shows the viscosity profiles at different time instants. Table 1. The flow-law parameters assumed for upper crust (plagioclase), lower crust (eclogite) and mantle (olivine). Mineral or rock  E (kJ mol−1)  n  A (MPa−ns−1)  References  Plagioclase  641  3  1012.7  Rybacki et al. (2006)  Eclogite  480  3.4  103.3  Jin et al. (2001)  Olivine  610  4.94  106.09  Korenaga & Karato (2008)  Mineral or rock  E (kJ mol−1)  n  A (MPa−ns−1)  References  Plagioclase  641  3  1012.7  Rybacki et al. (2006)  Eclogite  480  3.4  103.3  Jin et al. (2001)  Olivine  610  4.94  106.09  Korenaga & Karato (2008)  View Large The evolution of density can be approximated by the linear relation   $$\Delta \rho (z,t)=\alpha \rho _{0}(z)[T(z,0)-T(z,t)],$$ (48)where α is thermal expansivity and ρ0 is the reference potential density (Korenaga & Korenaga 2016, appendix A) defined at initial potential temperature (Tm). The thermal expansivity is assumed to be 3 × 10−5 K−1. The reference potential density is set to 3000, 3340, 3290 and 3320 kg m−3, respectively, for upper crust, lower crust, lithospheric mantle and asthenosphere (Korenaga 2006; Korenaga & Korenaga 2016; Johnson et al. 2014). Note that we should use potential density instead of in situ density. When a material is isentropically brought down to a greater depth, for example, its temperature and density both increase, so when comparing densities at different depths, we always need to take into account this density variation along the adiabat, which can be achieved using potential density. We estimate the evolution of such potential density with the cooling of the lithosphere (Fig. 3c). As the lower crust becomes negatively buoyant with respect to the underlying lithospheric mantle after ∼35 Myr, only the cases of 50 and 100 Myr exhibit the unstable density stratification in Fig. 3(c). We use the average density for each layer obtained from the density profile while calculating the dispersion relation. This approach minimized the number of primary density contrasts; we also computed the dispersion relation with continuous density variation with no significant difference. The density jump at the upper–lower crust interface is set to zero to take into account the effects of phase transition (Section 3.2). Utilizing the propagator matrix technique, the growth timescale is computed using the assumed viscosity and density profiles for a range of lithospheric ages (Fig. 3d). It can be seen that the growth timescale is initially very long (∼8 Gyr) when the eclogitic lower crust first appears at ∼35 Myr; this is because the dense layer is very thin at the beginning. As the layer thickness increases with cooling, the growth timescale drops to ∼200 Myr at 40 Myr, and then gradually increases to 10 Gyr after 100 Myr of cooling. This increase originates in the temperature dependence of viscosity. In addition to the continuous viscosity profile, we also use the minimum and average viscosities of each layer, the result of which forms an envelope covering all the likely values of the growth timescale for the reasonable range of lithospheric viscosity. The minimum timescale of such envelope turns out to be ∼150 Myr corresponding to a lithospheric age of 40–50 Myr; the growth time is simply too long to cause crustal delamination, because the timescale increases to a few billion years at the age of 100 Myr. In other words, when the eclogitic crust becomes thick enough to be gravitationally unstable, it is already too stiff to flow by cooling, as a consequence of viscosity being strongly dependent on temperature. In this particular model, therefore, the delamination is unlikely to occur by the Rayleigh–Taylor instability. 4.2 Growth time formula for general cases In order to study the delamination of the early crust in a comprehensive manner, we must explore the plausible ranges of the model parameters involved and devise a practical method to compute the growth timescale. The following parameter ranges are deemed sufficiently wide for our purpose: d1 = 5–50 km, d2 = 5–50 km, d3 = 50–200 km, Δρ2 = 10–100 kg m−3 and Δρ3 = −30 to 0 kg m−3. The age of the oceanic lithosphere (t) is varied from 1 to 100 Myr, and the range of the reference stress (σref) is chosen between 0.001 and 1 MPa. With a chosen reference stress, the viscosity profile may be constructed following the flow law (eq. 47) at each time instant, which can be used to compute the minimum growth timescale and corresponding wavelength. First, we consider a reference case with d1 = 25 km, d2 = 25 km, d3 = 120 km, Δρ2 = 50 kg m−3, Δρ3 = −15 kg m−3, t = 50 Myr and log (σref) = −2. This particular set of input parameters results a minimum growth timescale of ∼300 Myr at a perturbation wavelength of ∼100 km. Now the model space to be explored can be defined around this reference state with the following standard deviation:   $$\sigma (\mathbf {x}_{1})=[10,10,35,15,10,35,1]^{T},$$ (49)where   $$\mathbf {x}_{1}=[\Delta z_{1},\Delta z_{2},\Delta z_{3},\Delta \rho _{2},\Delta \rho _{3},t,\log (\sigma _{{\rm ref}})]^{T}.$$ (50)Here, Δzi and Δρi are the thickness of the ith layer and the density jump at the ith interface, respectively. This model space is subjected to random sampling using the normal distribution, and a model vector (eq. 50) is constructed for each sample. This obtained set of random model vectors is then used to compute the dispersion relation. An ensemble of minimum growth timescales is obtained through 1000 such calculations. Even though we use the continuously varying viscosity while calculating the growth rate, we summarize a viscosity profile by calculating the log-mean viscosity of each layer, and we use such log-mean viscosities in the following linear regression. That is, instead of using the age of the oceanic crust and reference stress as two model parameters, we use a set of log-mean viscosities, so that we can directly measure the sensitivity of growth rate to viscosity. The log-mean viscosities corresponding to the reference state are as follows: η1 = 1024 Pa s, η2 = 1023 Pa s, η3 = 1020 Pa s and η4 = 1019 Pa s, respectively. The new model vector is defined with the log-mean viscosities ηi:   $$\mathbf {x}_{2}=[\Delta z_{1},\Delta z_{2},\Delta z_{3},\Delta \rho _{2},\Delta \rho _{3},\log (\eta _{1}),\log (\eta _{2}),\log (\eta _{3}),\log (\eta _{4})]^{T},$$ (51)with the standard deviation of   \begin{eqnarray} \sigma (\mathbf {x}_{2})&=&[10,10,35,15,10,2,1,0.5,0.5]^{T}. \end{eqnarray} (52)Despite the complexity of our model, the minimum growth timescale and corresponding wavelength of perturbations are found to be predicted with reasonable accuracy using linear functions:   \begin{eqnarray} \log _{10}(\tau )&=&a_{0}+\sum _{i=1}^{9}a_{i}X_{i}, \end{eqnarray} (53)  \begin{eqnarray} \log _{10}(\lambda )&=&b_{0}+\sum _{i=1}^{9}b_{i}X_{i}, \end{eqnarray} (54)where Xi is the mean subtracted input parameter that is normalized by the corresponding standard deviation (eq. 52). Constants a0 through b9 are estimated by the least-squares regression: a0 = 2.4350, b0 = 2.0407, a1 = 0.1197, a2 = 0.6143, a3 = 0.0093, a4 = 0.0897, a5 = 0.0058, a6 = 0.1099, a7 = 0.8303, a8 = 0.1319, a9 = 0.00042, and b1 = 0.0027, b2 = 0.0304, b3 = 0.0013, b4 = 0.0024, b5 = 0.0037, b6 = 0.0344, b7 = 0.086, b8 = 0.2115 and b9 = 0.00026. Fig. 4 demonstrates a reasonable correspondence between the predicted and actual values of the model outputs for the aforementioned 1000 simulations. This way of summarizing simulation results enables us not only to see the sensitivity of the growth timescale on the model parameters, but also to quickly reproduce the important results without doing the whole calculation. Keeping only the model parameters with high enough sensitivities, we may suggest the following expressions for the growth timescale and the wavelength of perturbations:   $$\log _{10}(\tau )= 2.435-0.061(\Delta z_{2}-25)-0.009(\Delta \rho _{2}-50)+0.830\left(\log _{10}(\eta _{2})-23\right) + 0.264\left(\log _{10}(\eta _{3})-20\right),$$ (55)  \begin{eqnarray} \log _{10}(\lambda )&=&2.041+0.003(\Delta z_{2}-25)-0.034\left(\log _{10}(\eta _{1})-24\right) \end{eqnarray} (56)  \begin{eqnarray} &&\quad \quad\quad \quad-0.086\left(\log _{10}(\eta _{2})-23\right)-0.423\left(\log _{10}(\eta _{3})-20\right), \end{eqnarray} (57)where Δz, Δρ, η, τ and λ are in km, kg m−3, Pa s, Myr and km, respectively. Figure 4. View largeDownload slide Correspondence between the predicted and actual values for (a) minimum growth timescale and (b) the wavelength of perturbations. Figure 4. View largeDownload slide Correspondence between the predicted and actual values for (a) minimum growth timescale and (b) the wavelength of perturbations. While the model space is high-dimensional, the growth timescale is primarily a function of the eclogite layer thickness, the density jump at crust–mantle interface, and most importantly, the viscosity of the eclogite layer. These formulae for growth timescale and wavelength of perturbation are helpful to assess the likelihood of crustal delamination in the early Earth. Based on thermomechanical calculations, Johnson et al. (2014) suggested the possibility of crustal delamination by the Rayleigh–Taylor instability with the timescale of a few million years. Although they have treated the thermodynamics of the basalt–eclogite transition in a detailed manner, they seemed to have neglected the impact of crustal rheology on the instability. As we have shown, the growth timescale varies strongly with the temperature-dependent viscosity of eclogite; cold and stiff lower crust becomes resistant to flow, and the growth timescale may be too high to allow the possibility of delamination. Note that eclogite is weaker than pure garnet aggregates (Karato et al. 1995; Jin et al. 2001), but even with the rheology of eclogite, the lower crust is too strong. It is possible, however, to reduce the effective viscosity of eclogite by invoking higher reference stress. The likelihood of crustal delamination thus seems to be limited to special tectonic settings. 5 DISCUSSION The propagator matrix formulation of the Rayleigh–Taylor instability has been tested and validated against analytical solutions. In this semi-analytical technique, the instability problem is reduced to a compact eigenvalue problem, whose size depends only on the number of primary density jumps. Therefore, the method can handle a range of different situations, including continuous variation of material properties. Also because of the problem being reduced to finding out the eigenvalue of a reasonably small matrix, the computational complexity of the problem (i.e. O(M3)) is considerably lower than that of the existing methods (e.g. Bassi & Bonnin 1988; Conrad & Molnar 1997). This computational efficiency enables us to execute a large number of computations quickly. Thus, we can explore the high-dimensional model space extensively to delineate the sensitivity of the growth rate on various model parameters. In other words, this matrix-based technique allows us to gain efficiently useful insight into the physics of the multiple-layer Rayleigh–Taylor instability in general. Though linear stability analysis does not provide any information on finite-amplitude phenomena, we can sweep the whole parameter space quickly and identify the model subspace that deserves further investigation. The propagator matrix formulation together with the Monte Carlo sensitivity analysis may be applied to other geophysical situations as well. For example, massive terrestrial planets with stagnant lid convection may feature a temperature profile that enters into the stability field of the eclogite after crust grows beyond a critical thickness (O'Rourke & Korenaga 2012). Thus, we may expect gravitational instability at the crust–mantle interface. Accessing the likelihood of crustal delamination in such cases through our analysis may lead to a strong conclusion whether these planets can escape the stagnant lid regime. Also, eclogite may be produced through pressure-induced phase transition in thick continental crust (Sobolev & Babeyko 2005), making it negatively buoyant with respect to the underlying mantle. The delamination of such thick continental crust (e.g. Bird 1979; Meissner & Mooney 1998; Jull & Kelemen 2001) may also be examined using our approach. Our analysis of early crust delamination underscores the importance of rheology to the dynamics of delamination. Therefore, it is of paramount importance to better understand lithospheric viscosity if we aim to be more conclusive about the likelihood of delamination by gravitational instabilities. With further deformation experiments combined with proper statistical analysis (e.g. Korenaga & Karato 2008; Mullet et al. 2015), the rheology of crust and mantle materials is expected to be constrained with smaller uncertainties. The Monte Carlo sensitivity analysis can point to the subset of material properties that are most important from the geodynamical perspective, thereby facilitating the feedback between mineral physics and geodynamics. ACKNOWLEDGEMENTS This study is based on work supported by the NASA through the NASA Astrobiology Institute under cooperative agreement no. NNA15BB03A issued through the Science Mission Directorate. We thank Gael Choblet and Neil Ribe for their careful reviews and constructive comments. REFERENCES Bassi G., Bonnin J., 1988. Rheological modelling and deformation instability of lithosphere under extension, Geophys. J. Int.  93 485– 504. Google Scholar CrossRef Search ADS   Bird P., 1979. Continental delamination and the Colorado Plateau, J. geophys. Res.  84 7561– 7571. Google Scholar CrossRef Search ADS   Canright D., Morris S., 1993. Buoyant instability of a viscous film over a passive fluid, J. Fluid Mech.  255 349– 372. Google Scholar CrossRef Search ADS   Chandrasekhar S., 1961. Hydrodynamic and Hydromagnetic Stability , Clarendon, Oxford. Chu X., Korenaga J., 2012. Olivine rheology, shear stress, and grain growth in the lithospheric mantle: geological constraints from the Kaapvaal craton, Earth. planet. Sci. Lett.  333 52– 62. Google Scholar CrossRef Search ADS   Conrad C.P., Molnar P., 1997. The growth of Rayleigh-Taylor type instabilities in the lithosphere for various rheological and density structures, Geophys. J. Int.  129 95– 112. Google Scholar CrossRef Search ADS   Fletcher R.C., Hallet B., 1983. Unstable extension of the lithosphere: a mechanical model for Basin and Range structure, J. geophys. Res.  88 7457– 7466. Google Scholar CrossRef Search ADS   Foley S.F., Buhre S., Jacob D.E., 2003. Evolution of the Archean crust by delamination and shallow subduction, Nature  421 249– 252. Google Scholar PubMed  Forte A.M., 2000. Seismic-geodynamic constraints on mantle flow: implications for layered convection, mantle viscosity, and seismic anisotropy in the deep mantle, in Earth's Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale , eds Karato S., Forte A., Liebermann R., Masters G., Stixrude L., Geophys. Monogr. Vol. 117, pp. 3– 36. AGU, Washington, DC. Google Scholar CrossRef Search ADS   Hager B.H., O'Connell R.J., 1981. A simple global model of plate dynamics and mantle convection, J. geophys. Res.  86 4843– 4867. Google Scholar CrossRef Search ADS   Herzberg C., Condie K., Korenaga J., 2010. Thermal history of the Earth and its petrological expression, Earth planet. Sci. Lett.  292 79– 88. Google Scholar CrossRef Search ADS   Houseman G.A., Molnar P., 1997. Gravitational (Rayleigh-Taylor) instability of a layer with non-linear viscosity and convective thinning of continental lithosphere, Geophys. J. Int.  128 125– 150. Google Scholar CrossRef Search ADS   Ida S., Nakagawa Y., Nakazawa K., 1987. The Earth's core formation due to the Rayleigh-Taylor instability, Icarus  69 239– 248. Google Scholar CrossRef Search ADS   Jin Z.-M., Zhang J., Green H., Jin S., 2001. Eclogite rheology: implications for subducted lithosphere, Geology  29 667– 670. Google Scholar CrossRef Search ADS   Johnson T.E., Brown M., Kaus B.J., VanTongeren J.A., 2014. Delamination and recycling of Archean crust caused by gravitational instabilities, Nature Geosci.  7 47– 52. Google Scholar CrossRef Search ADS   Jull M., Kelemen P.B., 2001. On the conditions for lower crustal convective instability, J. geophys. Res.  106 6423– 6446. Google Scholar CrossRef Search ADS   Karato S.-i., 2012. Deformation of Earth Materials: An Introduction to the Rheology of Solid Earth  Reprint edn, Cambridge University Press. Karato S.-i., Wang Z., Liu B., Fujino K., 1995. Plastic deformation of garnets: systematics and implications for the rheology of the mantle transition zone, Earth. planet. Sci. Lett.  130 13– 30. Google Scholar CrossRef Search ADS   Korenaga J., 2006. Archean geodynamics and the thermal evolution of Earth, in Archean Geodynamics and Environments  eds Benn K., Mareschal J-C., Condie K.C., Geophys. Monogr. Vol. 164, pp. 7– 32. AGU, Washington, DC. Google Scholar CrossRef Search ADS   Korenaga J., Karato S.-I., 2008. A new analysis of experimental data on olivine rheology, J. geophys. Res.  113 B02403, doi:10.1029/2007JB005100. Google Scholar CrossRef Search ADS   Korenaga T., Korenaga J., 2016. Evolution of young oceanic lithosphere and the meaning of seafloor subsidence rate, J. geophys. Res.  121 6315– 6332. Google Scholar CrossRef Search ADS   Meissner R., Mooney W., 1998. Weakness of the lower continental crust: a condition for delamination, uplift, and escape, Tectonophysics  296 47– 60. Google Scholar CrossRef Search ADS   Mullet B.G., Korenaga J., Karato S.-I., 2015. Markov chain Monte Carlo inversion for the rheology of olivine single crystals, J. geophys. Res.  120 3142– 3172. Google Scholar CrossRef Search ADS   O'Rourke J.G., Korenaga J., 2012. Terrestrial planet evolution in the stagnant-lid regime: size effects and the formation of self-destabilizing crust, Icarus  221 1043– 1060. Google Scholar CrossRef Search ADS   Philpotts A., Ague J., 2009. Principles of Igneous and Metamorphic Petrology  Cambridge University Press. Google Scholar CrossRef Search ADS   Rybacki E., Gottschalk M., Wirth R., Dresen G., 2006. Influence of water fugacity and activation volume on the flow properties of fine-grained anorthite aggregates, J. geophys. Res.  111 B03203, doi:10.1029/2005JB003663. Google Scholar CrossRef Search ADS   Selig F., 1965. A theoretical prediction of salt dome patterns, Geophysics  30 633– 643. Google Scholar CrossRef Search ADS   Sobolev S.V., Babeyko A.Y., 2005. What drives orogeny in the Andes?, Geology  33 617– 620. Google Scholar CrossRef Search ADS   Stevenson D., 1990. Fluid dynamics of core formation, in Origin of the Earth , pp. 231– 249, eds Newsom H., Jones J.E., Oxford University Press, London. Taylor G., 1950. The instability of liquid surfaces when accelerated in a direction perpendicular to their planes, Proc. R. Soc. Lond. A  201 192– 196. Google Scholar CrossRef Search ADS   Turcotte D., Schubert G., 2002. Geodynamics  Cambridge University Press. Google Scholar CrossRef Search ADS   Whitehead J.A., Luther D.S., 1975. Dynamics of laboratory diapir and plume models, J. geophys. Res.  80 705– 717. Google Scholar CrossRef Search ADS   Zaleski S., Julien P., 1992. Numerical simulation of Rayleigh-Taylor instability for single and multiple salt diapirs, Tectonophysics  206 55– 69. Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations