Add Journal to My Library
Geophysical Journal International
, Volume 212 (3) – Mar 1, 2018

9 pages

/lp/ou_press/the-rayleigh-taylor-instability-in-a-self-gravitating-two-layer-e0sR3YSVjd

- Publisher
- The Royal Astronomical Society
- Copyright
- © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
- ISSN
- 0956-540X
- eISSN
- 1365-246X
- D.O.I.
- 10.1093/gji/ggx507
- Publisher site
- See Article on Publisher Site

Summary The dispersion relation of the Rayleigh–Taylor instability in the spherical geometry is of profound importance in the context of the Earth's core formation. Here we present a complete derivation of this dispersion relation for a self-gravitating two-layer viscous sphere. Such relation is, however, obtained through the solution of a complex transcendental equation, and it is difficult to gain physical insights directly from the transcendental equation itself. We thus also derive an empirical formula to compute the growth rate, by combining the Monte Carlo sampling of the relevant model parameter space with linear regression. Our analysis indicates that the growth rate of Rayleigh–Taylor instability is most sensitive to the viscosity of inner layer in a physical setting that is most relevant to the core formation. Core, Numerical solutions, Planetary interiors 1 INTRODUCTION The Rayleigh–Taylor instability of a self-gravitating fluid sphere is important in the context of planet formation, especially in Earth's core formation (e.g. Ida et al. 1989; Stevenson 1990; Hoogenboom & Houseman 2006; Sasaki & Abe 2007). A dynamical situation typically envisioned for an early, growing Earth can be characterized by an inner layer of highly viscous silicate materials and an outer layer of denser molten metal (e.g. Honda et al. 1993; Rubie et al. 2015). This gravitationally unstable structure is susceptible to the Rayleigh–Taylor instability at the interface of these two layers. Such instability may lead to the dripping of the metallic layer into the silicate layer, which is an important part of core formation. However, the growth timescale of this Rayleigh–Taylor instability depends on several parameters such as viscosity, density and layer thicknesses, and an accurate dispersion relation becomes desirable. The dispersion relation for the Rayleigh–Taylor instability in a two-layer self-gravitating sphere was published previously by Ida et al. (1989) for both inviscid and viscous cases. This dispersion relation was utilized to discuss the Earth's core formation (Ida et al. 1987). Whereas the inviscid case is correctly presented, however, their derivation for the viscous case is incorrect. In this research note, we derive a correct dispersion relation for the Rayleigh–Taylor instability in a two-layer self-gravitating viscous sphere. As the growth rate depends on a number of parameters, it can be challenging to understand the variation of the growth rate in the high-dimensional parameter space. Also, the estimation of the growth rate in case of finite Prandtl number is based on the solution of a non-linear transcendental equation, which can be very cumbersome in the spherical geometry. Therefore, we develop a practical approach to derive an approximate formula for the growth rate by combining Monte Carlo sampling and linear regression for l = 2, the mode most relevant to the Earth's core formation. The l = 1 mode is the translational mode and results simply in an overall shift of the centre of mass of the whole spherical system. Despite being the fast growing mode for the case of μ1 ≫ μ2, therefore, it is not relevant to the core formation (e.g. Elsasser 1963; Honda et al. 1993). Modes with l > 1 are important for core formation, but with increasing l, the viscous dissipation becomes higher in proportion to ≈l(l + 1), thereby making growth timescale longer. The l = 2 mode thus has the maximum contribution, and we restrict ourselves to this particular case when deriving an approximate formula. The structure of this note is following. We begin with the theoretical formulation of the Rayleigh–Taylor instability in a self-gravitating two-layer viscous sphere. We then present a few representative dispersion curves for different combinations of layer viscosities. We next discuss the practical approach to obtain an approximate formula for growth rate. Our result suggests that when the viscosity of the inner layer is much higher than that of the outer layer, the growth rate is controlled mostly by the properties of the inner layer as well as the density contrast of the two layers. 2 DERIVATION OF THE DISPERSION RELATION The governing equations for an incompressible viscous fluid include the equation of continuity \begin{equation} \nabla \cdot {\bf U}=0, \end{equation} (1)the equation of motion \begin{equation} \frac{{\rm d}\mathbf {U}}{{\rm d}t}=\nabla \cdot \mathbf {\Sigma }+{\bf F}, \end{equation} (2)and the constitutive relation \begin{equation} \mathbf {\Sigma }=-P\mathbf {I}+\mu \left[\nabla \mathbf {U}+(\nabla \mathbf {U})^{T}\right], \end{equation} (3)where U is the total velocity field, Σ is the total stress, F is the body force, P is the total pressure, I is the identity operator, 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 U}&=&{\bf U}_{0}+{\bf u}, \end{eqnarray} (4) \begin{eqnarray} \mathbf {\Sigma }&=&\mathbf {\Sigma }_{0}+\mathbf {\sigma }, \end{eqnarray} (5) \begin{eqnarray} P&=&P_{0}+\delta p, \end{eqnarray} (6) \begin{eqnarray} {\bf F}&=&{\bf F}_{0}+{\bf f}, \end{eqnarray} (7)where [U0, Σ0, P0, F0]T and [u, σ, δp, f]T represent the hydrostatic reference state (U0 is trivially zero) and the perturbed state, respectively. In the hydrostatic state, we may write the body force (gravity) as a gradient of potential that satisfies the Poisson equation: \begin{eqnarray} \mathbf {F}_{0}&=&-\rho _{0}\nabla \Psi , \end{eqnarray} (8) \begin{eqnarray} \nabla ^{2} \Psi &=& -4\pi G \rho _{0}, \end{eqnarray} (9)where ρ0 and G are the density of the reference state and the gravitational constant, respectively. We may solve the Poisson equation for the hydrostatic state (eq. 9) and obtain the unperturbed gravitational potential as \begin{eqnarray} \Psi &=& \frac{4\pi G}{3}\left(\rho _{1}(r^{2}/2-3R_{1}^{2}/2)+3\rho _{2}(R_{1}^{2}-R_{2}^{2})/2\right), \quad \quad \forall r<R_{1}, \end{eqnarray} (10) \begin{eqnarray} \quad&=&\frac{4\pi G}{3}\left(\rho _{1}(r^{2}/2-3R_{2}^{2}/2)+(\rho _{2}-\rho _{1})R_{1}^{3}/r\right), \quad \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (11)where ρ1, ρ2, R1 and R2 are the density of the inner layer, the density of the outer layer, the radius of the inner sphere, and the radius of the outer sphere, respectively (Fig. 1). The hydrostatic state confirms the equality of the pressure gradient and the gravitational force: \begin{equation} \nabla P_{0}=\rho _{0}\nabla \Psi , \end{equation} (12)which will be used in boundary conditions. Subtracting the hydrostatic reference state and neglecting higher order terms lead to \begin{eqnarray} \nabla _{j}u_{j}&=&0, \end{eqnarray} (13) \begin{eqnarray} \nabla _{j}\sigma _{ij}+f_{i}&=&\frac{\partial u_{i}}{\partial t}, \end{eqnarray} (14) \begin{eqnarray} \sigma _{ij}&=&-\delta p\delta _{ij}+\mu \left(\nabla _{j}u_{i}+\nabla _{i}u_{j}\right), \end{eqnarray} (15)where fi arises solely due to density perturbation at the interfaces of two adjacent fluids, and it can also be represented as the gradient of a potential perturbation δΨ satisfying the Poisson equation: \begin{eqnarray} f_{i}&=&-\rho _{0}\nabla \delta \Psi , \end{eqnarray} (16) \begin{eqnarray} \nabla ^{2} \delta \Psi &=&-4\pi G \delta \rho , \end{eqnarray} (17)where δρ is the perturbation to the reference density. Together with the perturbed eqs (13)–(17), we may write the equation of motion as \begin{equation} \frac{\partial \mathbf {u}}{\partial t}=-\nabla \Phi +\nu \nabla ^{2} \mathbf {u}, \end{equation} (18)where \begin{equation} \Phi =\frac{\delta P}{\rho }+\delta \Psi . \end{equation} (19)The system of eqs (13)–(19) describing the dynamics of the perturbed state is used to calculate the growth rate of initial density perturbations. In the 3-D spherical geometry, we expand the field variables in the spherical harmonics assuming the horizontal invariance of material properties (e.g. Chandrasekhar 1961) \begin{eqnarray} \Xi (r,\theta ,\phi ) &=& \sum _{l,m}\Xi _{lm}(r)Y^{m}_{l}(\theta ,\phi ). \end{eqnarray} (20)Using this separation of variables, we solve for the perturbed gravitational potential generated through the deformation of the interface (assuming the outer surface of the sphere to be rigid): \begin{eqnarray} \delta \Psi (r)&=&\left(\frac{4\pi G(\rho _{2}-\rho _{1})}{(2l+1)R_{1}^{l-1}}\delta r_{1}\right) r^{l}, \quad \quad \forall r < R_{1}, \end{eqnarray} (21) \begin{eqnarray} \qquad&=&\left(\frac{4\pi G (\rho _{2}-\rho _{1})}{2l+1}R_{1}^{l+2}\delta r_{1}\right)\bigg/r^{l+1}, \quad \quad \forall R_{1}<r<R_{2}. \end{eqnarray} (22)We can also readily solve for the total perturbed potential Φ by taking divergence of the eq. (18) and using the continuity eq. (13): \begin{equation} \nabla ^{2} \Phi =0, \end{equation} (23)the solution of which may further be used to solve the equation of motion (18). However, before proceeding to solve the equation of motion, it is convenient to perform poloidal–toroidal decomposition in the spherical geometry (e.g. Chandrasekhar 1961) as we will only require the poloidal part for linear stability analysis. Using such decomposition, we may express the perturbed velocity field as \begin{equation} \mathbf {u}=\sum _{l,m}\left(\nabla \times (T_{lm}(r)Y^{m}_{l}\hat{\mathbf {r}})+\nabla \times \left[\nabla \times (S_{lm}(r)Y^{m}_{l}\hat{\mathbf {r}})\right]\right), \end{equation} (24)where Slm and Tlm denote the poloidal and toroidal coefficients, respectively. Assuming u∝exp (σt) in the linear regime with σ being the growth rate, the poloidal part of the equation of motion (eq. 18) may be written as the following second-order differential equation: \begin{equation} \sigma S_{lm}=-\nabla \Phi +\nu \left[\frac{l(l+1)}{r^{2}}-D^{2}\right]S_{lm}, \end{equation} (25)where ν and D denote the kinematic viscosity and the radial derivative, that is, d/dr, respectively. The complementary function and particular solution of the eq. (25) can be obtained for the inner and outer layers separately as linear combinations of modified Bessel functions: \begin{eqnarray} S_{lm}(r)&=&A_{1}\sqrt{(}r)I_{l+1/2}(q_{1}r)-\frac{Q_{1}R_{1}^{2}}{\sigma }\left(\frac{r}{R_{1}}\right)^{l+1}, \quad \quad \forall r<R_{1},\nonumber \\ &=&A_{2}\sqrt{(}r)I_{l+1/2}(q_{2}r)+B_{2}\sqrt{r}K_{l+1/2}(q_{2}r)-\left[Q_{2}\left(\frac{r}{R_{1}}\right)^{l+1}+Q_{3}\left(\frac{r}{R_{1}}\right)^{-l}\right]R_{1}^{2}/\sigma , \quad \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (26)where q1 = σ/ν1, q2 = σ/ν2, ν1 and ν2 are the viscosities of the inner and outer layers, respectively, and [A1, Q1, A2, B2, Q2, Q3]T is a set of unknown constants to be obtained using the boundary conditions. The velocity field associated with the Rayleigh–Taylor instability may be obtained from the poloidal potential Slm. The radial (ur) and transverse (uθ) components are obtained as follows: \begin{eqnarray} u_{r}&=&\sum _{l,m}\frac{l(l+1)}{r^{2}}S_{lm}(r)Y_{l}^{m}(\theta ,\phi ),\nonumber \\ &=&\sum _{l,m}\left[\frac{A_{1}l(l+1)}{r^{3/2}}I_{l+1/2}(q_{1}r)-\frac{Q_{1}l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}\right]Y_{l}^{m}(\theta ,\phi ), \quad \forall r<R_{1}, \end{eqnarray} (27) \begin{eqnarray} &&{\quad=\sum _{l,m}\Bigg[\frac{A_{2}l(l+1)}{r^{3/2}}I_{l+1/2}(q_{2}r)+\frac{B_{2}l(l+1)}{r^{3/2}}K_{l+1/2}(q_{2}r)}\nonumber \\ &&{\quad\quad-\frac{l(l+1)}{r^{2}}\left[Q_{2}\left(\frac{r}{R_{1}}\right)^{l+1}+Q_{3}\left(\frac{r}{R_{1}}\right)^{-l}\right]R_{1}^{2}/\sigma \Bigg]Y_{l}^{m}(\theta ,\phi ), \quad \forall R_{1}<r<R_{2},} \end{eqnarray} (28) \begin{eqnarray} u_{\theta }&=&\sum _{l,m}\frac{1}{r}\frac{{\rm d}S_{lm}(r)}{{\rm d}r}\frac{\partial Y_{l}^{m}(\theta ,\phi )}{\partial \theta },\nonumber \\ &=&\sum _{l,m}\Big[A_{1}\left((l+1)r^{-3/2}I_{l+1/2}(q_{1}r)+q_{1}r^{-1/2}I_{l+3/2}(q_{1}r)\right)-\frac{Q_{1}(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}\Big]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall r<R_{1}, \end{eqnarray} (29) \begin{eqnarray} &&{\quad=\sum _{l,m}\Bigg[A_{2}\left((l+1)r^{-3/2}I_{l+1/2}(q_{2}r)+q_{2}r^{-1/2}I_{l+3/2}(q_{2}r)\right)+B_{2}\left((l+1)r^{-3/2}K_{l+1/2}(q_{2}r)-q_{2}r^{-1/2}K_{l+3/2}(q_{2}r)\right)}\nonumber \\ &&{\quad\quad-\frac{1}{\sigma }\left(Q_{2}(l+1)\left(\frac{r}{R_{1}}\right)^{l-1}-Q_{3}l \left(\frac{r}{R_{2}}\right)^{-(l+2)}\right)\Bigg]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall R_{1}<r<R_{2}.} \end{eqnarray} (30)The stresses may be calculated from the velocities using the following constitutive relations: \begin{eqnarray} \sigma _{rr}&=&2\mu \frac{\partial u_{r}}{\partial r}, \end{eqnarray} (31) \begin{eqnarray} \sigma _{r\theta }&=&\mu \left[\frac{1}{r}\frac{\partial u_{r}}{\partial \theta }-\frac{u_{\theta }}{r}+\frac{\partial u_{\theta }}{\partial r}\right]. \end{eqnarray} (32) Figure 1. View largeDownload slide The spherical system considered in the Rayleigh–Taylor instability analysis. Unstable stratification corresponds to ρ2 > ρ1. Figure 1. View largeDownload slide The spherical system considered in the Rayleigh–Taylor instability analysis. Unstable stratification corresponds to ρ2 > ρ1. It is essential to compute the stresses correctly, particularly the normal stress, as it is used in the discontinuity condition, which is directly related to the expression of the growth rate. The previously published result by Ida et al. (1989) uses an incorrect expression for the normal and tangential stresses (their eqs A-21 and A-22); a simple dimensional analysis reveals incompatibility between the both sides of these equations. To make calculations simpler, we may express the velocities and stresses in the following manner: \begin{eqnarray} u_{r}&=&\sum _{l,m}\left[a_{1}A_{1}+b_{1}Q_{1}\right]Y_l^m(\theta ,\phi ), \quad \forall r<R_{1}, \end{eqnarray} (33) \begin{eqnarray} \quad&=&\sum _{l,m}\left[a_{2}A_{2}+b_{2}B_{2}+c_{2}Q_{2}+d_{2}D_{2}\right]Y_l^m(\theta ,\phi ), \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (34) \begin{eqnarray} u_{\theta }&=&\sum _{l,m}\left[e_{1}A_{1}+f_{1}Q_{1}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall r<R_{1}, \end{eqnarray} (35) \begin{eqnarray} \quad=\sum _{l,m}\left[e_{2}A_{2}+f_{2}B_{2}+g_{2}Q_{2}+h_{2}Q_{3}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (36) \begin{eqnarray} \tau _{rr}&=&-P+2\mu _{1}\sum _{l,m}\left[i_{1}A_{1}+j_{1}Q_{1}\right]Y_l^m(\theta ,\phi ), \quad \forall r<R_{1}, \end{eqnarray} (37) \begin{eqnarray} \quad&=&-P+2\mu _{2}\sum _{l,m}\left[i_{2}A_{2}+j_{2}B_{2}+k_{2}Q_{2}+l_{2}Q_{3}\right]Y_l^m(\theta ,\phi ), \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (38) \begin{eqnarray} \tau _{r\theta }&=&\mu _{1}\sum _{l,m}\left[\alpha _{1}A_{1}+\beta _{1}Q_{1}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall r<R_{1}, \end{eqnarray} (39) \begin{eqnarray} \quad\quad&=&\mu _{2}\sum _{l,m}\left[\alpha _{2}A_{2}+\beta _{2}B_{2}+\gamma _{2}Q_{2}+\delta _{2}Q_{3}\right]\frac{\partial Y_{l}^{m}}{\partial \theta }, \quad \forall R_{1}<r<R_{2}, \end{eqnarray} (40)where \begin{eqnarray} a_{1}(r)&=&\frac{l(l+1)}{r^{3/2}}I_{l+1/2}(q_{1}r), \end{eqnarray} (41) \begin{eqnarray} b_{1}(r)&=&-\frac{l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}, \end{eqnarray} (42) \begin{eqnarray} a_{2}(r)&=&\frac{l(l+1)}{r^{3/2}}I_{l+1/2}(q_{2}r), \end{eqnarray} (43) \begin{eqnarray} b_{2}(r)&=&\frac{l(l+1)}{r^{3/2}}K_{l+1/2}(q_{2}r), \end{eqnarray} (44) \begin{eqnarray} c_{2}(r)&=&-\frac{l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}, \end{eqnarray} (45) \begin{eqnarray} d_{2}(r)&=&-\frac{l(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{-(l+2)}, \end{eqnarray} (46) \begin{eqnarray} e_{1}(r)&=&\frac{(l+1)}{r^{3/2}}I_{l+1/2}(q_{1}r)+\frac{q_{1}}{\sqrt{r}}I_{l+3/2}(q_{1}r), \end{eqnarray} (47) \begin{eqnarray} f_{1}(r)&=&-\frac{(l+1)}{\sigma }\left(\frac{r}{R_{1}}\right)^{l-1}, \end{eqnarray} (48) \begin{eqnarray} e_{2}(r)&=&(l+1)r^{-3/2}I_{l+1/2}(q_{2}r)+q_{2}r^{-1/2}I_{l+3/2}(q_{2}r), \end{eqnarray} (49) \begin{eqnarray} f_{2}(r)&=&(l+1)r^{-3/2}K_{l+1/2}(q_{2}r)-q_{2}r^{-1/2}K_{l+3/2}(q_{2}r), \end{eqnarray} (50) \begin{eqnarray} g_{2}(r)&=&-\frac{(l+1)}{\sigma }(r/R_{1})^{l-1}, \end{eqnarray} (51) \begin{eqnarray} h_{2}(r)&=&\frac{l}{\sigma }\left(\frac{r}{R_{2}}\right)^{-(l+2)}, \end{eqnarray} (52) \begin{eqnarray} i_{1}(r)&=&l(l+1)\left[(l-1)r^{-5/2}I_{l+1/2}(q_{1}r)+q_{1}r^{-3/2}I_{l+3/2}(q_{1}r)\right], \end{eqnarray} (53) \begin{eqnarray} j_{1}(r)&=&-\frac{l(l+1)(l-1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (54) \begin{eqnarray} i_{2}(r)&=&l(l+1)\left[(l-1)r^{-5/2}I_{l+1/2}(q_{2}r)+q_{2}r^{-3/2}I_{l+3/2}(q_{2}r)\right], \end{eqnarray} (55) \begin{eqnarray} j_{2}(r)&=&l(l+1)\left[(l-1)r^{-5/2}K_{l+1/2}(q_{2}r)-q_{2}r^{-3/2}K_{l+3/2}(q_{2}r)\right], \end{eqnarray} (56) \begin{eqnarray} k_{2}(r)&=&-\frac{l(l+1)(l-1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (57) \begin{eqnarray} l_{2}(r)&=&\frac{l(l+1)(l+2)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{-(l+3)}, \end{eqnarray} (58) \begin{eqnarray} \alpha _{1}(r)&=&2(l+1)(l-1)r^{-5/2}I_{l+1/2}(q_{1}r)+(2l+1)q_{1}r^{-3/2}I_{l+3/2}(q_{1}r)+q_{1}^{2}r^{-1/2}I_{l+5/2}(q_{1}r), \end{eqnarray} (59) \begin{eqnarray} \beta _{1}(r)&=&-\frac{2(l-1)(l+1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (60) \begin{eqnarray} \alpha _{2}(r)&=&2(l+1)(l-1)r^{-5/2}I_{l+1/2}(q_{2}r)+(2l+1)q_{2}r^{-3/2}I_{l+3/2}(q_{2}r)+q_{2}^{2}r^{-1/2}I_{l+5/2}(q_{2}r), \end{eqnarray} (61) \begin{eqnarray} \beta _{2}(r)&=&2(l+1)(l-1)r^{-5/2}K_{l+1/2}(q_{2}r)-(2l+1)q_{2}r^{-3/2}K_{l+3/2}(q_{2}r)+q_{2}^{2}r^{-1/2}K_{l+5/2}(q_{2}r), \end{eqnarray} (62) \begin{eqnarray} \gamma _{2}(r)&=&-\frac{2(l-1)(l+1)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{l-2}, \end{eqnarray} (63) \begin{eqnarray} \delta _{2}(r)&=&-\frac{2l(l+2)}{\sigma R_{1}}\left(\frac{r}{R_{1}}\right)^{-(l+3)}. \end{eqnarray} (64) We are now ready to use boundary conditions to obtain the unknown coefficients. We assume the outer boundary to be rigid, that is, ur = uθ = 0 at r = R2. The inner boundary is assumed to be welded and hence the field variables are continuous across the deformed boundary. The first-order Taylor expansion leads to the continuity of the velocity and transverse stress across the r = R1 boundary, that is, $$\left[u_{r}\right]_{R_{1}-}^{R_{1}+}=0$$, $$\left[u_{\theta }\right]_{R_{1}-}^{R_{1}+}=0$$, and $$\left[\sigma _{r\theta }\right]_{R_{1}-}^{R_{1}+}=0$$, while the radial stress is continuous only across the deformed boundary, $$\left[-(P_{0}+\delta P)+2\mu \partial _{r}u_{r}\right]_{R_{1}-\delta r}^{R_{1}+\delta r}=0$$, and becomes discontinuous across r = R1 due to the density jump. P0 and δP may be calculated from eqs (12), (16), and (23). The rigid and continuous boundary conditions together with the orthonormality of the spherical harmonics form the following complete set of linear equations of the unknown coefficients for each (l, m) pair: \begin{eqnarray} a_{2}(R_{2})A_{2}+b_{2}(R_{2})B_{2}+c_{2}(R_{2})Q_{2}+d_{2}(R_{2})Q_{3}&=&0, \end{eqnarray} (65) \begin{eqnarray} e_{2}(R_{2})A_{2}+f_{2}(R_{2})B_{2}+g_{2}(R_{2})Q_{2}+h_{2}(R_{2})Q_{3}&=&0, \end{eqnarray} (66) \begin{eqnarray} a_{1}(R_{1})A_{1}+b_{1}(R_{1})Q_{1}&=&a_{2}(R_{1})A_{2}+b_{2}(R_{1})B_{2}+c_{2}(R_{1})Q_{2}+d_{2}(R_{1})Q_{3}, \end{eqnarray} (67) \begin{eqnarray} e_{1}(R_{1})A_{1}+f_{1}(R_{1})Q_{1}&=&e_{2}(R_{1})A_{2}+f_{2}(R_{1})B_{2}+g_{2}(R_{1})Q_{2}+h_{2}(R_{1})Q_{3}, \end{eqnarray} (68) \begin{eqnarray} \mu _{1}\left[\alpha _{1}(R_{1})A_{1}+\beta _{1}(R_{1})Q_{1}\right]&=&\mu _{2}[\alpha _{2}(R_{1})A_{2}+\beta _{2}(R_{1})B_{2}+\gamma _{2}(R_{1})Q_{2}+\delta _{2}(R_{1})Q_{3}], \end{eqnarray} (69) \begin{eqnarray} \frac{4\pi G(\rho _{1}-\rho _{2})[2(l-1)\rho _{1}+3\rho _{2}]}{3(2l+1)}\delta r&=&(l+1)(\rho _{1}Q_{1}-\rho _{2}Q_{2})+l\rho _{2}Q_{3}+2\mu _{2}/R_{1}\partial _{r}u_{r2}-2\mu _{1}/R_{1}\partial _{r}u_{r1}. \end{eqnarray} (70)We may eliminate the several unknowns from these equations and obtain the following form: \begin{eqnarray} A_{2}&=&\xi _{1}Q_{2}+\eta _{1}Q_{3}, \end{eqnarray} (71) \begin{eqnarray} B_{2}&=&\xi _{2}Q_{2}+\eta _{2}Q_{3}, \end{eqnarray} (72) \begin{eqnarray} A_{1}&=&X_{1}Q_{2}+Y_{1}Q_{3}, \end{eqnarray} (73) \begin{eqnarray} Q_{1}&=&X_{2}Q_{2}+Y_{2}Q_{3}, \end{eqnarray} (74) \begin{equation} Q_{2}=\Theta Q_{3}, \end{equation} (75) \begin{equation} Q_{3}=\Sigma \delta r, \end{equation} (76)where \begin{eqnarray} \xi _{1}&=&-\frac{f_{2}(R_{2})c_{2}(R_{2})-b_{2}(R_{2})g_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (77) \begin{eqnarray} \eta _{1}&=&-\frac{f_{2}(R_{2})d_{2}(R_{2})-b_{2}(R_{2})h_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (78) \begin{eqnarray} \xi _{2}&=&-\frac{a_{2}(R_{2})g_{2}(R_{2})-e_{2}(R_{2})c_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (79) \begin{eqnarray} \eta _{2}&=&-\frac{a_{2}(R_{2})h_{2}(R_{2})-e_{2}(R_{2})d_{2}(R_{2})}{a_{2}(R_{2})f_{2}(R_{2})-e_{2}(R_{2})b_{2}(R_{2})}, \end{eqnarray} (80) \begin{eqnarray} X_{1}&=&\frac{f_{1}(R_{1})[a_{2}(R_{1})\xi _{1}+b_{2}(R_{1})\xi _{2}+c_{2}(R_{1})]-b_{1}(R_{1})[e_{2}(R_{1})\xi _{1}+f_{2}(R_{1})\xi _{2}+g_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (81) \begin{eqnarray} Y_{1}&=&\frac{f_{1}(R_{1})[a_{2}(R_{1})\eta _{1}+b_{2}(R_{1})\eta _{2}+d_{2}(R_{1})]-b_{1}(R_{1})[e_{2}(R_{1})\eta _{1}+f_{2}(R_{1})\eta _{2}+h_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (82) \begin{eqnarray} X_{2}&=&\frac{a_{1}(R_{1})[e_{2}(R_{1})\xi _{1}+f_{2}(R_{1})\xi _{2}+g_{2}(R_{1})]-e_{1}(R_{1})[a_{2}(R_{1})\xi _{1}+b_{2}(R_{1})\xi _{2}+c_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (83) \begin{eqnarray} Y_{2}&=&\frac{a_{1}(R_{1})[e_{2}(R_{1})\eta _{1}+f_{2}(R_{1})\eta _{2}+h_{2}(R_{1})]-e_{1}(R_{1})[a_{2}(R_{1})\eta _{1}+b_{2}(R_{1})\eta _{2}+d_{2}(R_{1})]}{a_{1}(R_{1})f_{1}(R_{1})-b_{1}(R_{1})c_{1}(R_{1})}, \end{eqnarray} (84) \begin{equation} \Theta =\frac{\mu _{2}[\alpha _{2}(R_{1})\eta _{1}+\beta _{2}(R_{1})\eta _{2}+\delta _{2}(R_{1})]-\mu _{1}[\alpha _{1}(R_{1})Y_{1}+\beta _{1}(R_{1})Y_{2}]}{\mu _{1}[\alpha _{1}(R_{1})X_{1}+\beta _{1}(R_{1})X_{2}]-\mu _{2}[\alpha _{2}(R_{1})\xi _{1}+\beta _{2}(R_{1})\xi _{2}+\gamma _{2}(R_{1})]}, \end{equation} (85) \begin{equation} \Sigma =\frac{4\pi G (\rho _{1}-\rho _{2})(2(l-1)\rho _{1}+3\rho _{2})}{3(2l+1)\Pi }, \end{equation} (86)and \begin{eqnarray} \Pi &=&\left[\left((l+1)\rho _{1}X_{2}-(l+1)\rho _{2}\right)+\frac{2\mu _{2}}{R_{1}}\left(i_{2}(R_{1})\xi _{1}+j_{2}(R_{1})\xi _{2}+k_{2}(R_{1})\right)-\frac{2\mu _{1}}{R_{1}}\left(i_{1}(R_{1})X_{1}+j_{1}(R_{1})X_{2}\right)\right]\Theta \nonumber\\ &&+\left[(l+1)\rho _{1}Y_{2}+l\rho _{2}+\frac{2\mu _{2}}{R_{1}}\left(i_{2}(R_{1})\eta _{1}+j_{2}(R_{1})\eta _{2}+l_{2}(R_{1})\right)-\frac{2\mu _{1}}{R_{1}}\left(i_{1}(R_{1})Y_{1}+j_{1}(R_{1})Y_{2}\right)\right]. \end{eqnarray} (87)We may write the relation between the amplitude of perturbation and the vertical velocity at the interface as \begin{eqnarray} \frac{{\rm d}(\delta r)}{{\rm d}t}&=&u_{r}(R_{1}). \end{eqnarray} (88)By noting that δr ∝ exp (σt) in the linear regime, eqs (34) and (88) together with the six unknowns (i.e. A1-Q3) may be combined to obtain the expression for the growth rate \begin{equation} \sigma =\left(a_{1}(R_{1})[X_{1}\Theta +Y_{1}]+b_{1}(R_{1})[X_{2}\Theta +Y_{2}]\right)\Sigma . \end{equation} (89)This is a transcendental equation for σ, because its right-hand side depends on σ in a complex manner. 3 EXAMPLES AND APPROXIMATE FORMULA In this section, we provide examples of dispersion relation for several combinations of layer viscosities. It is convenient to present the growth rate in the dimensionless form, so we scale the growth rate with respect to the free fall rate as follows: \begin{equation} \sigma^{\prime}=\sigma /\sqrt{4\pi G (\rho _{2}-\rho _{1})/3}. \end{equation} (90)First, we fix the viscosity of the inner layer (i.e. μ1) and compute the growth rate for a range of outer layer viscosity (i.e. μ2). Similarly, we compute the growth rates for a range of inner layer viscosity keeping the viscosity of the outer layer constant. In both of these cases, the primary contributing factor turns out to be max (μ1, μ2) (Fig. 2). For μ1 ≫ μ2, the translational mode approximately corresponds to the motion of a rigid sphere through a viscous fluid. The corresponding growth rate turns out to be the fastest; the inner rigid sphere almost instantaneously translates to a new position leading to an overall shift of the centre of mass. Viscosity variation therefore does not affect considerably the growth rate for this translational mode as long as μ1 ≫ μ2 (Figs 2a and b). The growth rates of modes with l ≥ 2 exhibit nearly linear dependence on the inner layer viscosity (Fig. 2b). In the case of μ2 ≫ μ1, little variation in growth rate is observed with μ1, while growth rates associated with modes with small angular order (i.e. l < 5) vary linearly with μ2 (Figs 2c and d). However, the case of μ2 ≫ μ1 is not relevant in the context of core formation. Figure 2. View largeDownload slide The dispersion relations for several combinations of inner and outer layer viscosities: (a) the case where μ2 is varied keeping μ1 constant and μ1 ≫ μ2, (b) the case where μ1 is varied keeping μ2 constant and μ1 ≫ μ2—most relevant to the core formation, (c) the case where μ1 is varied keeping μ2 constant and μ2 ≫ μ1, (d) the case where μ2 is varied keeping μ1 constant and μ2 ≫ μ1. Other parameters are fixed as follows: R1 = 106 m, dR = (R2 − R1) = 104.84 m, ρ1 = 103.74 kg m−3 and δρ = 103.5 kg m−3. Figure 2. View largeDownload slide The dispersion relations for several combinations of inner and outer layer viscosities: (a) the case where μ2 is varied keeping μ1 constant and μ1 ≫ μ2, (b) the case where μ1 is varied keeping μ2 constant and μ1 ≫ μ2—most relevant to the core formation, (c) the case where μ1 is varied keeping μ2 constant and μ2 ≫ μ1, (d) the case where μ2 is varied keeping μ1 constant and μ2 ≫ μ1. Other parameters are fixed as follows: R1 = 106 m, dR = (R2 − R1) = 104.84 m, ρ1 = 103.74 kg m−3 and δρ = 103.5 kg m−3. In order to study the Rayleigh–Taylor instability associated with core formation in a comprehensive manner, it would be convenient to derive an approximate formula for the growth rate for the most important mode associated with the deformation of the interface. To derive such a formula, we explore the following parameter ranges: R1 = 105 to 107 m, ΔR = (R2 − R1) = 103 to 5 × 106 m, ρ1 = 3 × 103 to 104 kg m−3, δρ = (ρ2 − ρ1) = 103 to 104 kg m−3, μ1 = 1015 to 1025 Pa s, and μ2 = 106 to 1010 Pa s. First, we consider a reference case with R1 = 106 m, ΔR = 104.84 m, ρ1 = 103.74 kg m−3, δρ = 103.5 kg m−3, μ1 = 1020 Pa s, and μ2 = 108 Pa s. The model space to be explored can be defined around this reference case with the following standard deviation: \begin{equation} \zeta (\mathbf {x})=[0.5, 0.5,0.5,1,3,1.5]^{T}, \end{equation} (91)where \begin{equation} \mathbf {x}=[\log _{10}(R_{1}),\log _{10}(\rho _{1}),\log _{10}(\delta \rho ),\log _{10}(\Delta R),\log _{10}(\mu _{1}),\log _{10}(\mu _{2})]^{T}. \end{equation} (92)This model space is subjected to random sampling using the normal distribution, and a model vector (eq. 92) is constructed for each sample. This obtained set of random model vectors is then used to compute the growth rate corresponding to the l = 2 mode. An ensemble of growth rates is obtained through 200 such calculations. Despite the complexity of the transcendental equation (eq. 89), the growth rate is found to be predicted with reasonable accuracy using the following function: \begin{equation} \log _{10}(\sigma )=a_{0}+\sum _{j=1}^{6}a_{j}X_{j}, \end{equation} (93)where Xj is the mean subtracted input parameter that is normalized by the corresponding standard deviation (eq. 91). Constants a0 through a6 are estimated by least square regression: a0 = −10.43, a1 = 1.86, a2 = 0.87, a3 = 1.3, a4 = 0.047, a5 = −0.9752, and a6 = −0.054. Fig. 3 demonstrates a reasonable correspondence between the predicted and actual values of growth rates for the aforementioned 200 computations. This way of summarizing enables us not only to see the sensitivity of the growth rate on the model parameters, but also to quickly compute the growth rate without doing the whole calculation. We suggest the following expression for the growth rate corresponding to the l = 2 mode: \begin{equation} \sigma =10^{-10.13}R_{1}^{1.86}(\rho _{2}-\rho _{1})^{1.3}\rho _{1}^{0.87}(R_{2}-R_{1})^{0.047}\mu _{1}^{-0.98}\mu _{2}^{-0.054}, \end{equation} (94)where R, ρ, μ and σ are in m, kg m−3, Pa s and s−1, respectively. Figure 3. View largeDownload slide Correspondence between the actual and predicted values of growth rate for l = 2. Figure 3. View largeDownload slide Correspondence between the actual and predicted values of growth rate for l = 2. 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 Chandrasekhar S., 1961. Hydrodynamic and Hydromagnetic Stability Clarendon. Elsasser W., 1963. Early history of the Earth, in Earth Science and Meteoritics pp. 1– 30, eds Geiss J., Goldberg E., North-Holland. Honda R., Mizutani H., Yamamoto T., 1993. Numerical simulation of Earth's core formation, J. geophys. Res. 98 2075– 2089. Google Scholar CrossRef Search ADS Hoogenboom T., Houseman G.A., 2006. Rayleigh–Taylor instability as a mechanism for corona formation on Venus, Icarus 180 292– 307. 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 Ida S., Nakagawa Y., Nakazawa K., 1989. The Rayleigh–Taylor instability in a self-gravitating two-layer fluid sphere, Earth Moon Planets 44 149– 174. Google Scholar CrossRef Search ADS Rubie D., Nimmo F., Melosh H., 2015. Formation of the Earth's core, in Treatise On Geophysics 2nd edn, vol. 9, pp. 51– 90, ed. Stevenson D., Elsevier. Sasaki T., Abe Y., 2007. Rayleigh–Taylor instability after giant impacts: imperfect equilibration of the Hf-W system and its effect on the core formation age, Earth Planets Space 59 1035– 1045. 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 Univ. Press. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Geophysical Journal International – Oxford University Press

**Published: ** Mar 1, 2018

Loading...

personal research library

It’s your single place to instantly

**discover** and **read** the research

that matters to you.

Enjoy **affordable access** to

over 12 million articles from more than

**10,000 peer-reviewed journals**.

All for just $49/month

Read as many articles as you need. **Full articles** with original layout, charts and figures. Read **online**, from anywhere.

Keep up with your field with **Personalized Recommendations** and **Follow Journals** to get automatic updates.

It’s easy to organize your research with our built-in **tools**.

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.

## “Hi guys, I cannot tell you how much I love this resource. Incredible. I really believe you've hit the nail on the head with this site in regards to solving the research-purchase issue.”

Daniel C.

## “Whoa! It’s like Spotify but for academic articles.”

@Phil_Robichaud

## “I must say, @deepdyve is a fabulous solution to the independent researcher's problem of #access to #information.”

@deepthiw

## “My last article couldn't be possible without the platform @deepdyve that makes journal papers cheaper.”

@JoseServera

- Read unlimited articles
- Personalized recommendations
- No expiration
- Print 20 pages per month
- 20% off on PDF purchases
- Organize your research
- Get updates on your journals and topic searches

**$49/month**

Start Free Trial

14-day Free Trial

Best Deal — 39% off
### Annual Plan

- All the features of the Professional Plan, but for
**39% off**! - Billed annually
- No expiration
- For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

~~$588~~

**$360/year**

Start Free Trial

14-day Free Trial

Read and print from thousands of top scholarly journals.

System error. Please try again!

or

By signing up, you agree to DeepDyve’s Terms of Service and Privacy Policy.

Already have an account? Log in

Bookmark this article. You can see your Bookmarks on your DeepDyve Library.

To save an article, **log in** first, or **sign up** for a DeepDyve account if you don’t already have one.