3-D numerical modelling of crustal polydiapirs with volume-of-fluid methods

SUMMARY

Gravitational instabilities exert a crucial role on the Earth dynamics and in particular on its differentiation. The Earth’s crust can be considered as a multilayered fluid with different densities and viscosities, which may become unstable in particular with variations in temperature. With the specific aim to quantify crustal scale polydiapiric instabilities, we test here two codes, JADIM and OpenFOAM, which use a volume-of-fluid (VOF) method without interface reconstruction, and compare them with the geodynamics community code ASPECT, which uses a tracking algorithm based on compositional fields. The VOF method is well-known to preserve strongly deforming interfaces. Both JADIM and OpenFOAM are first tested against documented two and three-layer Rayleigh–Taylor instability configurations in 2-D and 3-D. 2-D and 3-D results show diapiric growth rates that fit the analytical theory and are found to be slightly more accurate than those obtained with ASPECT. We subsequently compare the results from VOF simulations with previously published Rayleigh–Bénard analogue and numerical experiments. We show that the VOF method is a robust method adapted to the study of diapirism and convection in the Earth’s crust, although it is not computationally as fast as ASPECT. OpenFOAM is found to run faster than, and conserve mass as well as JADIM. Finally, we provide a preliminary application to the polydiapiric dynamics of the orogenic crust of Naxos Island (Greece) at about 16 Myr, and propose a two-stages scenario of convection and diapirism. The timing and dimensions of the modelled gravitational instabilities not only corroborate previous estimates of timing and dimensions associated to the dynamics of this hot crustal domain, but also bring preliminary insight on its rheological and tectonic contexts.

Geomechanics, Numerical modelling, Crustal structure, Diapirism, Dynamics: gravity and tectonics, Rheology: crust and lithosphere

1 INTRODUCTION

Thermally and chemically driven gravitational instabilities are the main processes involved in the differentiation of the Earth (Schubert & Turcotte 1971; Christensen 1984; Christensen & Yuen 1985; Korenaga 2018). In turn, the Earth’s crust itself is subject to differentiation owing to magmatism, metamorphism and deformation (Taylor & McLennan 1985; Rudnick & Fountain 1995). It has been proposed that partial melting of orogenic roots modifies their rheology and allows for the development of gravitational instabilities that play a key role in controlling crustal differentiation (Ramberg 1980, 1981a; Perchuk et al. 1992; Burg & Vanderhaeghe 1993; Brown 1994; Sawyer 1994; Weinberg & Podladchikov 1994; Cruden et al. 1995; Jull & Kelemen 2001; Vigneresse 2006; Gerya et al. 2008; Vanderhaeghe 2009; Gerbault et al. 2018).

Gneiss domes observed worldwide are interpreted as the result of gravitational instabilities developed within a hot partially molten crust (Whitney et al. 2004). Active convection of partially molten crust has also been proposed, involving several tens of cubic kilometres of weak and light material during generally more than 10 Myr (Riel et al. 2016; Vanderhaeghe et al. 2018). As an example case, the area of Naxos island (Greece) presents typical domes and subdomes which have been interpreted as possibly resulting from polydiapirism and convection (Vanderhaeghe et al. 2018).

Whereas several studies have simulated gravitational instabilities of two or more crustal layers such as Poliakov et al. (1993) or Wilcock & Whitehead (1991), very few have, to our knowledge, precisely compared their results with theory. Ramberg (1972) compared his theory with geophysical data, while Berner et al. (1972) compared numerical and experimental models with Ramberg’s theory Ramberg (1981a), explaining their misfit by an insufficient number of elements and an ‘oversized’ time step. Some studies have predicted semi-analytically the growth rate and wavelength for two layer systems: in an infinite half-space as Whitehead (1988), Kaus (2004), Burg et al. (2004), Schmalholz & Podladchikov (1999) or in a finite space (Selig 1965; Turcotte & Schubert 1982). However as a matter of fact, polydiapiric crustal structures do not appear to have been much studied independently from a context of tensile or compressional tectonic drive, since the analytical and experimental models of Rayleigh–Taylor instabilities carried out by Ramberg (1981a) and the numerical models by Weinberg & Schmeling (1992). Field observations of polydiapirism leave open the question of whether relatively small dome structures first form independently at depth and then merge or coalesce to form a single larger structure (dome) close to the surface, or whether this large dome already develops at depth and then progressively develops additional instabilities upon rising (progressive segregation and subdomes). Finding out the correct process would allow to identify from which depth and which environment do specific elements and compositions appear, separate or co-exist, constraining the evolution of distinct elementary compositions, some of which lead to the concentration of mineral resources (Toé et al. 2013; Eglinger et al. 2016; Menant et al. 2018).

On the other hand, convection in the Earth’s mantle has been extensively studied both numerically and experimentally [Bercovici & Schubert (2009) and references there in]. In this well developed field, numerical codes can encounter issues with mass conservation and numerical diffusion at layers interfaces (Deubelbeiss & Kaus 2008; Schmeling et al. 2008; Hillebrand et al. 2014; Heister et al. 2017; Pusok et al. 2017). This is an even more critical issue when dealing with crustal scale instabilities, since crustal differentiation involves the recurrent separation, coalescence and segregation of layers of highly contrasted compositions that evolve during melting and deformation. These segregation and coalescence processes occur at different scales, and require a robust numerical tool in terms of tracking the evolution of chemical interfaces at the intermediate scale of a few hundred metres. This motivates the use of a volume-of-fluid (VOF) method, dedicated to the conservation of chemical interfaces. Recently, Puckett et al. (2018) compared a VOF method (using interface reconstruction) against other standard methods in models of mantle convection. They qualitatively studied two-layer Rayleigh–Taylor and Rayleigh–Bénard systems in 2-D, and concluded that their VOF method may be the most appropriate method for modelling interfaces separating chemical compositions. In the present contribution, we present two existing codes namely JADIM and OpenFOAM built on VOF methods without interface reconstruction, and which were initially developed for other purposes (bubbles, drops, free-surface flows). JADIM is an in-house Fortran code developed at IMFT, and was proved to accurately describe two- and three-layer flows involving strongly deforming interfaces (Bonometti & Magnaudet 2006; Bonhomme et al. 2012). OpenFOAM in turn, is an open-source C++ code that offers a VOF solver for fluid mechanics (Jasak et al. 2007). OpenFOAM has been used in the Geosciences community in the last years (Orgogozo et al. 2014; Dietterich et al. 2017) but with other solvers adapted to the modelling of water and lava flows. We here test both codes for the development of crustal scale convective instabilities in 2-D and 3-D.

We aim at checking the accuracy and performance of both JADIM and OpenFOAM codes when modelling Rayleigh–Taylor and Rayleigh–Bénard instabilities in two and three dimensions. We therefore compare the computed solutions with available analytical solutions as well as with those obtained with the open source mantle convection code ASPECT version 2.1.0 (Kronbichler et al. 2012; Heister et al. 2017; Bangerth et al. 2019). We first check whether the VOF method without interface reconstruction is able to reproduce complex 2-D and 3-D Rayleigh–Taylor and Rayleigh–Bénard crustal scale systems. Based on the accuracy and performance of the results, one code is chosen to model gravitational instabilities in the specific context of Naxos Island (Greece). As such the present contribution stands as a preliminary study, that first validates the VOF method when used to model crustal flows, and second paves the way to a second contribution aiming at further exploring the influence of subscale crustal properties, then fully benefiting from the robustness of the VOF method.

This study is structured as follows. We first introduce the numerical methods (Section 2). The results for standard iso-thermal two- and three-layer systems (Rayleigh–Taylor instability) are then presented in Section 3, compared with reference models (Weinberg & Schmeling 1992; van Keken et al. 1997) and linear stability theory. 3-D simulations are then compared with the 2-D three-layer systems. In Section 4, Rayleigh–Bénard instabilities in one- and two-layer systems are simulated and compared with experimental results (Le Bars & Davaille 2004; Vatteville et al. 2009). In Section 5, we discuss the codes performances and their suitability to model polydiapirism. Section 6 presents a preliminary application to Naxos, in which we illustrate the macroscopic thermomechanical setting with which the Naxos crust would have been able to develop its characteristic domes and sub-domes (Vanderhaeghe et al. 2018).

2 NUMERICAL METHODS

2.1 The VOF method

The VOF method is a fixed-grid approach based on the one-fluid model and considers that the various immiscible fluids (or ‘phases’) can be described as a single fluid whose local physical properties, namely density and viscosity, vary in space and time depending on the volume fraction C_i of each phase i (Hirt & Nichols 1981; Youngs 1982). The volume fraction of each fluid intrinsically obeys |$\sum \limits _{{i=1}}^n C_i = 1$| where n is the number of phases. In this study, we consider one-, two- and three-phase systems and 1 ≤ n ≤ 3. Typically, C_i = 1 in grid cells filled only with fluid i, and 0 < C_i < 1 in grid cells cross-cut by an interface. At this point need is to mention the two main classes of VOF methods: methods that try to reconstruct exactly the interface between fluids (e.g. Puckett et al. 2018), which requires significant computational time, and methods that do not, such as in the present approaches with JADIM and OpenFOAM. With no interface reconstruction, the thickness of the interfacial region is defined by 0 < C_i < 1, and typically occupies two to three grid cells. Dimensionless equations are presented in Appendix  B. Without any thermal effect, the local density ρ and viscosity μ of the fluids follow the relations:

$$\begin{eqnarray*} \rho = \sum \limits _{{i=1}}^n C_i \rho _i; \quad \mu = \left(\sum \limits _{{i=1}}^n \frac{C_i}{\mu _i} \right)^{-1}. \end{eqnarray*}$$

(1)

Note that the first equation in (1) is exact while the second equation is an ad hoc approximation [details given further below for the choice of μ’s interpolation, and see Schmeling et al. (2008)]. When thermal effects are taken into account (Section 4), the properties of the fluid are expressed as a function of both the volume fraction C_i and the temperature T, as specified below (eqs 6 and 7).

The general set of solved governing equations are the transport of the volume fractions, mass, momentum and energy equations, expressed as:

$$\begin{eqnarray*} \frac{\partial C_i}{\partial t } + {\bf U} \cdot \nabla C_i = - \nabla \cdot ({\bf U}_r C_r ), \end{eqnarray*}$$

(2)

$$\begin{eqnarray*} \nabla \cdot {\bf U} = 0, \end{eqnarray*}$$

(3)

$$\begin{eqnarray*} \rho \frac{\partial {\bf U}}{\partial t } + \rho {\bf U} \cdot \nabla {\bf U} = - \nabla P + \skew{2}\tilde{\rho }{\bf g} + \nabla \cdot \left[ \mu (\nabla {\bf U} + (\nabla {\bf U})^T)\right], \end{eqnarray*}$$

(4)

$$\begin{eqnarray*} \frac{\partial T}{\partial t} + {\bf U} \cdot \nabla T = \nabla \cdot (\kappa \nabla T), \end{eqnarray*}$$

(5)

where U, P, T are the velocity, pressure and temperature of the flow respectively, and |$\textbf {g}$| is the gravitational acceleration. It is worth noting that when used, the temperature-dependency of density is only applied to the gravitational term |$\tilde{\rho }{\bf g}$| (Boussinesq approximation). In particular, we set

$$\begin{eqnarray*} \skew{2}\tilde{\rho }= \sum \limits _{{i=1}}^n C_i \rho _i \mathcal {F}_i (T) ; \quad \tilde{\mu }= \left(\sum \limits _{{i=1}}^n \frac{C_i}{ \mu _i \mathcal {G}_i (T)} \right)^{-1}, \end{eqnarray*}$$

(6)

$$\begin{eqnarray*} \mathcal {F}_i (T) = 1 - \alpha _i (T-T^{ref}_i), \end{eqnarray*}$$

(7)

where |$T^{ref}_i$| and α_i are the reference temperature and the thermal expansion coefficient of fluid i, respectively, and |$\mathcal {G}_i (T)$| is function of temperature. Since |$\mathcal {G}_i (T)$| is model-specific, it will be defined below in each case.

In theory the r.h.s. of (2) should be zero, which is the case for JADIM and ASPECT. However, for OpenFOAM the term −∇ · (U_rC_r) is artificially added to (2) to reduce the effects of numerical smearing of the interface, where C_r = C₁ · (1 − C₁), and U_r, designated by Berberović et al. (2009) as a ‘compression velocity’, is evaluated at cell faces as a volume flux based on the maximum velocity magnitude in the interface region. This velocity is obtained from a face interpolation using the normalized variable diagram of Jasak et al. (1999, NVD approach), who proposed a high resolution differentiation scheme with a limiter (based on the ratio between volumetric flux gradients calculated at adjacent cell faces and cell centres). We refer the interested reader to Berberović et al. (2009) for further information.

Eqs (2)–(5) are solved on a structured staggered grid (JADIM) or on a complex mixture of collocated/staggered grid (OpenFOAM) using a finite-volume technique and a projection method used to enforce incompressibility. Domain decomposition and MPI parallelization is performed to allow for high-resolution simulations with a large number of grid cells.

2.2 Solvers characteristics

In JADIM, eq. (2) is solved using a modified version of the Flux Corrected Transport technique of Zalesak (1979) while eqs (4)–(5) are solved via a third order Runge–Kutta/Crank–Nicolson time-advancement scheme, the spatial gradients being approximated by a second-order central difference discretization. Pressure is computed by solving a Poisson pseudo-equation using the PETSC library. The linear system is solved by a Jacobi preconditioned conjugate gradient technique. The full numerical approach has been extensively described by Bonometti & Magnaudet (2007) and is not repeated here. JADIM is available only in the frame of a collaboration with IMFT.

In OpenFOAM (version 4.0 is used here), eq. (2) is solved using a Multidimensional Universal Limiter with Explicit Solution (MULES, Deshpande et al. (2012)) while eqs (3)–(4) are solved by a Pressure Implicit Splitting Operator (PISO) algorithm (Issa 1986), in which the momentum equation is solved first, considering the pressure at the previous time step. The solution gives an approximation of the new velocity field, which is used to solve the pressure equation to provide a first estimate of the new pressure field. This loop is then repeated with this new pressure field in order to correct the new velocity field, until a pre-determined tolerance is reached. Eq. (5) is solved by a first order Euler implicit time stepping algorithm. The OpenFOAM solvers named InterFoam, Multiphase-InterFoam (Rayleigh–Taylor instability) are used for Section 3, and are directly available in OpenFOAM. To solve the Rayleigh–Bénard cases of Section 4, we built our own solver which combines InterFoam and BuoyantBoussinesq-PimpleFoam [a two-fluid VOF solver of the energy eq. (5)], available here: https://gitlab.com/AurelieLN/openfoam.git.

ASPECT (version 2.1.0 in optimized mode is used here, Kronbichler et al. 2012; Heister et al. 2017; Bangerth et al. 2019) is a versatile state-of-the-art Finite Element code which is based on the deal.II library (Arndt et al. 2019). It is currently one of the most used open-source codes developed by the geodynamics community studying mantle convection and/or lithospheric deformation. Taylor-Hood elements are used for the Stokes equations (Donea & Huerta 2003), while second order elements are used for the temperature equation. The numerical methods implemented in ASPECT relevant to this study are presented in Kronbichler et al. (2012) and Heister et al. (2017). Material tracking is based on compositional fields (one field per fluid) and their corresponding advection equations are solved with second-order elements too, stabilized with the Entropy Viscosity method (Guermond et al. 2011). Although ASPECT allows the user (i) to tune the parameters of the Entropy Viscosity method, (ii) to use the SUPG method (Brooks & Hughes 1982) instead, (iii) to use a discontinuous Galerkin Finite Element method (He et al. 2017) for the advection–diffusion equations and (iv) to use the Particle-in-Cell method (Gassmöller et al. 2018, 2019), we felt that such explorations are outside the scope of this work and therefore resorted to using only default parameters. While a VOF method has been implemented in ASPECT (Robey & Puckett 2019), we do not use it here, first because it was published subsequently to the main development of our study, and second because at the time of writing it is only available in 2-D. Table 1 summarizes the main features of JADIM, OpenFOAM and ASPECT.

Table 1.

Numerical characteristics of the codes used in this work.

	JADIM	OpenFOAM	ASPECT
General method	VOF	VOF	Compositional fields
Grid	Eulerian	Eulerian	Eulerian
Time	3rd order Runge-Kutta/Crank-Nicolson	1st order Euler implicit	BDF-2
Spatial gradient	2nd order central difference	upwind	Q2 + EV stabilization

2.3 Inertia and Archimedes number

JADIM and OpenFOAM solve the Navier–Stokes equation on a Cartesian grid, and their formulation gives the possibility to take into account inertial processes (without extra cost in computational performance). In the considered Rayleigh–Taylor problems (Weinberg & Schmeling 1992; van Keken et al. 1997) inertia is negligible in eq. (4), and its presence in the formulation does not affect the results. Indeed, eqs (128) and (131) of Chandrasekhar (2013) and Ramberg (1968) showed that for a two-layer Rayleigh–Taylor system with a single kinematic viscosity ν, inertia can be neglected when the wavelength |$\displaystyle {\lambda \leqslant 2 \pi \left(\frac{4\nu ^2}{At \, g} \right)^{1/3}}$|⁠, with |$\displaystyle {At = \frac{\rho _1-\rho _2}{\rho _1+\rho _2}}$|⁠, and this is the case for the systems studied here.

Furthermore, since both JADIM and OpenFOAM are not fully implicit in the viscosity term, our numerical time steps are constrained to remain small. To circumvent unrealistic computational times, we decrease the viscosity artificially and show that this does not modify the characteristics of the modelled instabilities. To show this, we use the Archimedes number which compares buoyancy forces to viscous forces:

$$\begin{eqnarray*} \displaystyle {Ar = \frac{ \rho _1 (\rho _1-\rho _2) g H^3 }{\mu _1^2}} \end{eqnarray*}$$

(8)

We show in Appendix  C that a system with Ar <<1 and a system with Ar = 1 display equivalent behaviours. This behaviour can be assimilated to the independence of the corresponding drag coefficient on the Reynolds number in the Stokes flow regime (Duan et al. 2015).

3 MULTILAYER RAYLEIGH–TAYLOR INSTABILITIES

We define a series of numerical experiments in 2-D and 3-D Cartesian geometries to test the accuracy and the efficiency of the VOF method. In this section, the density and viscosity of each phase are kept constant and independent of temperature. Thus only eqs (2)–(4) are solved with |$\mathcal {F}(T)=\mathcal {G}(T)=1$|⁠. We first show the standard van Keken et al. (1997) two-layer system and compare the results produced by JADIM, OpenFOAM and ASPECT. Then, we study a three-layer system using the setup defined by Weinberg & Schmeling (1992) for crustal scale polydiapirism. The original results are compared with those computed with JADIM, OpenFOAM and ASPECT. We display a comparison with linear stability theory in all cases and show the evolution of the mass error. Below, we first briefly describe the linear stability theory and then present results for two-layer and three-layer systems in 2-D. The last subsection deals with 3-D three-layer systems.

3.1 Linear stability theory (LST)

The linear stability theory of Ramberg (1981a) predicts the dispersion relation between the wavelength and the growth rate of interfaces for a system with two or more horizontal layers of fluid of arbitrary density and viscosity, in an unstable configuration. The position y_i of the interface that delimits each layer i takes the form:

$$\begin{eqnarray*} y_i(t) = A_i e^{\kappa _i q_i t}, \end{eqnarray*}$$

(9)

where A_i depends on the initial perturbation, σ_i = κ_iq_i is the growth rate and q_i is defined as

$$\begin{eqnarray*} \displaystyle {q_i = \frac{(\rho _i - \rho _{i+1})g h_{i+1}}{2 \mu _{i+1}}}. \end{eqnarray*}$$

(10)

The parameter κ_i is a cumbersome combination of the dimensionless parameters |$\displaystyle {\gamma = \mu _i/\mu _{i+1}}$|⁠, |$\displaystyle {h_i/H}$|⁠, |$\displaystyle {\rho _i/\rho _{i+1}}$| and |$\displaystyle {(\rho _i-\rho _{i+1})(\rho _{i+1}-\rho _{i+2})}$|⁠, and its detailed formulation is found by solving the linear system of Ramberg (1981a, see eq. A1). Note that the theory presented in Ramberg (1981a) is valid when the interface height does not move more than 10 per cent of the wavelength.

Time, length and velocity are scaled by |$\mathcal {T}$|⁠, |$\mathcal {L}$| and |$\mathcal {U}$| as follows, using layer 2 as reference

$$\begin{eqnarray*} \mathcal {T}=q^{-1} = \frac{2 \mu _2}{(\rho _1-\rho _2) g h_2}; \quad \mathcal {L}=h_2; \quad \mathcal {U}=qh_2. \end{eqnarray*}$$

(11)

3.2 Two-layer Rayleigh–Taylor system

3.2.1 Numerical setup

We first start with the standard test case of a two-layer viscous Rayleigh–Taylor instability proposed by van Keken et al. (1997). A fluid of density ρ₁ and viscosity μ₁ is located above a less dense fluid of density ρ₂, viscosity μ₂ and thickness h₂ in a rectangular box of length L and height H. The two fluids are separated by a sinusoidal interface of low amplitude (Fig. 1). No-slip boundary conditions are used along the horizontal walls, free-slip conditions are used along the vertical walls, and a zero normal gradient is imposed for the volume fraction. For our numerical models considering a two layer system, we set Ar = 1, ρ₁/ρ₂ = 1.1, h₂/H = 0.2, and vary γ = μ₁/μ₂ = [1, 10, 10²].

Figure 1.

(a) Initial geometry of the two-layer Rayleigh–Taylor system and boundary conditions used in the van Keken et al. (1997) test case: h₂/H = 0.2 and L/H = 0.9142. (b) Initial geometry of the three-layer Rayleigh–Taylor system and boundary conditions used by Weinberg & Schmeling (1992): γ = μ₁/μ₂ and μ₂ = μ₃, h₂ = h₃ (see text for the specific value of the height ratios and length to height ratios).

In van Keken et al. (1997) a comparison between various numerical approaches was made (e.g. finite-element vs. finite-difference methods using either tracers, markers or field functions for the interface motion). We choose as reference results, those obtained with the finite-element method and a marker chain named Pvk. The reference grid has a resolution of 91 × 100 elements along the horizontal and vertical direction, respectively. The time evolution of several quantities is monitored, in particular:

• The dimensionless rms velocity: |$\displaystyle {V_{rms}/\mathcal {U} = \frac{1}{\mathcal {U}}\sqrt{\frac{1}{\mathcal {A}} \int _{\mathcal {A}} \Vert {\bf U}\Vert ^2 dS } }$|⁠, where |$\mathcal {A}=L\times H$| is the area of the computational domain.

• The dimensionless growth rate for each layer, cf. eq. (10): |$\displaystyle {\kappa = \frac{\sigma }{q}}$| (the slope of the V_rms curve at the beginning of the system destabilization),

• The relative velocity error: |$\displaystyle {\Delta V_{rms} = \frac{V_{rms} -V_{rms} ^{ref}}{V_{rms}^{ref}}}$|⁠, where |$V^{ref}_{rms}$| is a reference velocity defined later,

• The relative mass error of a phase M with respect to its initial mass M₀: |$\displaystyle {\Delta M = \frac{M-M_0}{M_0}}$|⁠. ΔM is the integral of one fluid within its own area versus the area it occupied at t = 0.

Averaging of the viscosity over an interface is known to potentially lead to different results (Schmeling et al. 2008; Bangerth et al. 2019). Most common methods are the arithmetic averaging (⁠|$\mu = \sum \limits _{{i=1}}^n C_i { \mu _i}$|⁠), the geometric averaging (⁠|$\mu = \mu _1^{C_1} \mu _2^{C_2}$|⁠), and the harmonic averaging (⁠|$\mu =\left(\sum \limits _{{i=1}}^n \frac{C_i}{ \mu _i} \right)^{-1}$|⁠). In van Keken et al. (1997), the method for viscosity averaging over interfaces is not specified. However in a following publication (Tackley & King 2003), this comparison is used and commented with an arithmetic averaging. Moreover, if we compare our results obtained using the arithmetic averaging (Figs D1 andD2) and those using the harmonic averaging ( Figs   D2 andD3), we observe that the arithmetic averaging gives results closer to those of van Keken et al. (1997). In fact, for a light viscous sphere rising in a more viscous fluid, the arithmetic averaging is found to be more appropriate than the harmonic averaging (Benkenida 1999). In the context of this benchmark the situation is alike, with a less viscous blob rising within a more viscous fluid.

3.2.2 Comparison of geometric patterns

We display in Fig. 2 snapshots of the time evolution of the two-phase system for three codes, JADIM, OpenFOAM and ASPECT for the viscosity ratio γ = 100 which displays the most drastic differences (other viscosity ratios are presented in appendix, Fig. D1). We superpose two types of visualization: isocontours of phase fractions of fluid 2 and colours. Black zones correspond to a phase fraction of fluid 2 between 0.66 and 1, grey zones to a phase fraction between 0.33 and 0.66, and white zones to a phase fraction between 0 and 0.33. The grey zones (where the interface diffuses) are broader for ASPECT than for the two VOF codes. It is only a visual and qualitative comparison.

Figure 2.

Time evolution of the two-phase system for a viscosity ratio of γ = 100, arithmetic viscosity averaging, at |$t/\mathcal {T}=50, 100, 150$| (columns).Van Keken's snapshots are modified from Fig. 6 of van Keken et al. (1997).The grid size is 91 × 100 for JADIM, OpenFOAM, and ASPECT. Iso-contours of the volume fraction of fluid 2 (initial bottom layer) are C₂ = 0.05, 0.5, 0.95 for JADIM, OpenFOAM and ASPECT. White zone: 0 < C₂ < 0.33; grey zone: 0.33 < C₂ < 0.66; black zone: 0.66 < C₂ < 1.

When the light fluid (fluid 2, black layer) reaches the top right corner of the domain, a layer of heavy fluid remains attached to the top wall (fluid 1, white layer), and continues to drip slowly. The general location of this layer’s interface is similar in all three codes to that of van Keken et al. (1997). However, looking in detail at this layer, we see that all codes develop distinct interface shapes. For instance with OpenFOAM and ASPECT, the secondary central drip plunges faster already at t/τ = 50 than with JADIM or van Keken et al. (1997).

For the test cases γ = 1 and γ = 10 in turn, the location of interfaces remains more similar for all three codes (see Appendix, Fig. D1). Note that the results obtained with OpenFOAM display small instabilities at the bottom of the model domain (Fig. D1b), similar to what Tackley & King (2003) had obtained with relatively high resolution and a great number of tracers. This is not the case for JADIM since the numerical thickness of the interface is slightly larger (two-three grid cells) than in OpenFOAM. At time t/τ = 150, a bubble is rising with OpenFOAM, is forming with JADIM, but does not exist in van Keken et al. (1997).

In conclusion, it is delicate at this stage to determine which method should produce the most correct results, and we can only state that structures formed with ASPECT display more diffusion than JADIM or OpenFOAM, which will be corroborated by the following numerical sensitivity tests.

3.2.3 Numerical sensitivity of the results

We plot in Fig. 3a) the time evolution of the dimensionless root mean square velocity (⁠|$V_{rms}/\mathcal {U}$|⁠), for different viscosity ratios. For γ = 1, the first velocity peak (inset A in Fig. 3a) is always well described even at low resolution (22 × 25) for both JADIM and OpenFOAM. van Keken et al. (1997) had observed that the arrival of the second diapir somewhat depends on the numerical code, and in fact the second diapir produced by OpenFOAM rises slightly faster than that simulated by van Keken et al. (1997) and JADIM (by about t/τ = 10), but it displays the same maximum value inset B in (Fig. 3a).

Figure 3.

(a and b) Comparison of the time evolution of the mean velocity in the Rayleigh–Taylor two-layer system by van Keken et al. (1997), with JADIM and OpenFOAM (arithmetic viscosity averaging and resolution 91 × 100). Velocity is scaled by |$\mathcal {U}$| and time is scaled by |$\mathcal {T}$| (eq. 11), for various viscosity contrasts |$\displaystyle {\gamma = \mu _1/\mu _2}$|⁠. A and B zoom on different scales at subsequent times. (b) Effect of the spatial resolution on velocity V_rms (scaled by |$\mathcal {U}$|⁠) at t/τ = 0.0981 for JADIM, OpenFOAM and ASPECT (γ = 1). N_y is the number of cells along the y-axis. The general trend is a line of slope ≈2.

The largest discrepancy occurs for γ = 100 where the results of JADIM and OpenFOAM underestimate (by 18 and 14 per cent, respectively) the ascent velocity predicted by van Keken et al. (1997) (Fig. 3a). For a comparison, we also plot in Fig. 3(a) the linear stability solution (Ramberg 1981a, see Appendix A1 for more details). For all viscosity ratios, simulations growth rates differ by at most 4 per cent from the linear stability theory.

Fig. 3(c) presents the relative error of all codes with respect to their own reference value against N_y, the number of grid cell along the y-axis. The reference value of each code is the mean velocity V_rms obtained at t/τ = 0.0981 with the finest grid. First, relative error decreases as the number of cells increases, showing grid convergence. Secondly, for a given resolution, the relative error is similar for codes JADIM and OpenFOAM (for example, for 1/N_y = 0.02, |$\Delta V_{\rm{rms}} = 0.05\,{\rm{per\,cent}} $|⁠) and generally lower for ASPECT (for 1/Ny = 0.02, |$\Delta V_{\rm{rms}} = 0.004\,{\rm{per\,cent}}$|⁠). Third, the slope of the curve is roughly 2, indicating that these approaches are second-order accurate in space.

The time evolution of the relative mass error ΔM (in per cent) for JADIM, OpenFOAM and ASPECT is presented in Fig. G1. For JADIM and OpenFOAM, ΔM oscillates but is always below 10⁻⁶ which is very small. ΔM for ASPECT is significantly larger compared to the VOF codes, with ΔM ≈ 10⁻².

These results show that the two VOF methods JADIM and OpenFOAM reproduce well the van Keken et al. (1997) test case. Tackley & King (2003) had previously explained differences in between models as resulting from the viscosity interpolation scheme with time stepping and with mesh resolution. We show here in addition, that VOF methods such as JADIM and OpenFOAM produce less interface diffusion and less mass error than ASPECT’s composition method.

3.3 Three-layer Rayleigh–Taylor systems

We now consider Rayleigh–Taylor instabilities in three-layer systems, which have been previously studied by Ramberg (1981a) and Weinberg & Schmeling (1992). The problem setup is similar to that of the previous section with the addition of a third fluid layer (see Fig. 1b). Weinberg & Schmeling (1992) studied various configurations of density and viscosity varying within these layers, assuming that each may represent a different composition or contain variable proportions of molten crustal material. They used a finite-difference method to solve the Stokes equations and tracked the interfaces using markers. In what follows, we study three of the configurations proposed by Weinberg & Schmeling (1992), and compare the results produced by JADIM, OpenFOAM, and for one case by ASPECT (case II), with the linear stability analysis of Ramberg (1981a) (briefly recalled in Section 3.1 for multilayers). We extrapolate these setups to 3-D in the following section.

3.3.1 RT numerical setups: three cases

Following Weinberg & Schmeling (1992), three configurations are tested. In case I, only the bottom layer 3 is gravitationally unstable. In cases II and III, both the intermediate and bottom layers 2 and 3 are gravitationally unstable. Since we showed in Section 2.3 and Appendix  C that the dynamics of a system remain similar when Ar ≤ 1, our simulations have Ar ∼ 1 while in Weinberg & Schmeling (1992) Ar ∼ 10⁻³³. More specifically, our physical parameters are identical to those of Weinberg & Schmeling (1992) except for the viscosities which are scaled smaller. The length of the computational domain is set to L/H = 2.24 for case I and L/H = 2.4 for cases II and III, so that it can include at least one dominant wavelength (Weinberg & Schmeling 1992). Both interfaces are perturbed by a random perturbation. In the present computations, our grid size is 224 × 100 for case I and 240 × 100 for cases II and III along the horizontal and vertical direction, respectively.

The general three-layer problem is described by five parameters, namely two Archimedes numbers based on properties of fluid 1 and 3, the viscosity ratio, the density ratio and the height ratio. We set μ₂ = μ₃ and h₂ = h₃ and density ratios close to one (they vary in each case), as in Weinberg & Schmeling (1992). Therefore, the problem reduces to the dimensionless parameters defined in section 3.2, namely, Ar ≈ 1, γ = μ₁/μ₂ = μ₁/μ₃ = [100, 200] and h₂/H = h₃/H = 0.125.

3.3.2 RT numerical results: comparison of geometrical patterns

Fig. 4 displays the time evolution of the three-layer Rayleigh–Taylor instability for cases I, II and III. Note that time is scaled by |$\mathcal {T}$| (eq. 11) and in order to make a relevant comparison, the origin of time t = 0 has been arbitrarily set as the time when one of the interfaces rises by a minimum height of 3 × 10⁻⁵H (temporal offset values are given for each 2-D and 3-D simulations in the figure captions). Indeed, we insert an initial perturbation at the interfaces (as in Weinberg & Schmeling 1992) which may not be defined exactly in the same manner for each numerical code depending on how the interface is calculated. For JADIM and OpenFOAM, interfaces at heights h₂ and h₃ are defined with a random perturbation of 0.24 per cent the local volume fraction of the cell crossed by the interface [as in Weinberg & Schmeling (1992)]. For ASPECT, both interfaces oscillate as z ≥ h_i + (rand() − 0.5) × 2.4 10⁻¹⁰). Thus, the moment of destabilization occurs at different times for each code.

Figure 4.

Three-phase Rayleigh–Taylor system by Weinberg & Schmeling (1992). Weinberg and Schmeling's snapshots are modified from Figs 3 and 5 of Weinberg & Schmeling (1992). Time evolution of the layers for the different codes: (a) Case I, the middle layer is the densest and γ = 100, (b) Case II, the top layer is the densest and γ = 100, and (c) Case III, the top layer is the densest and γ = 200. Time is scaled by |$\mathcal {T}$| with origin |$t/\mathcal {T}=0$| set at the time when one of the interfaces raised by a distance of 3 × 10⁻⁵H, Time offsets are, for case I: |$\text{Off}_{JA}^1 = 26$|⁠, |$\text{Off}_{OP}^1 = 24$|⁠, case II: |$\text{Off}_{JA}^2 = 300$|⁠, |$\text{Off}_{OP}^2 = 276$| and case III: |$\text{Off}_{JA}^3 = 564$|⁠, |$\text{Off}_{OP}^3 = 552$|⁠.

(a) Case I: ρ₂ > ρ₁ > ρ₃, γ = 100,Ar = 0.5

Here, ρ₁ = 0.89ρ₂ = 1.18ρ₃. The interface separating fluid 2 and fluid 3, from now on denoted as ‘interface 2-3’, deforms and forms regularly spaced rising domes. These domes merge inside layer 2 to form bigger domes which then rise through layer 1. Weinberg & Schmeling (1992)’s simulation displays 7 diapirs against 8.5 for JADIM and 6.5 for OpenFOAM.

(b) Case II: ρ₁ > ρ₂ > ρ₃, γ = 100, Ar = 0.981

Here, ρ₁ = 1.11ρ₂ = 1.12ρ₃. A mushroom shaped diapir composed of fluid 3 surrounded by fluid 2 (intricated domes structure) rises through layer 1. Diapirs rising to the left for JADIM and OpenFOAM are very similar to the diapir of Weinberg & Schmeling (1992) although it rises at the centre of the modelled domain. The diapirs rising to the right for JADIM and OpenFOAM present smaller wavelengths of the interface 2–3 than the diapir of Weinberg & Schmeling (1992). Nevertheless, the overall shape of the diapirs simulated by JADIM and OpenFOAM is very similar.

As for ASPECT, only one diapir forms at the centre of the domain, as in Weinberg & Schmeling (1992), however the interface 2–3 develops more diapirs (around five at t/τ = 441), which differs from Weinberg & Schmeling (1992) but is rather similar to those obtained with JADIM and OpenFOAM. In fact, the results of ASPECT are very sensitive to the random perturbation inserted at the interfaces and which depends directly on mesh resolution. ASPECT results resemble more those of Weinberg & Schmeling (1992) when the mesh is refined in the vertical direction (150 points with respect to 100, see Fig. E2c).

(c) Case III: ρ₁ > ρ₂ > ρ₃, γ = 200, Ar = 0.9

Here, ρ₁ = 1.1ρ₂ = 1.12ρ₃. The viscosity contrast between the uppermost and the intermediate layers is larger in this case. As a consequence the instability of interface 1–2 develops more slowly than in the previous cases, and the intermediate fluid sinks down rather than rise. Results obtained with JADIM and OpenFOAM are similar but the agreement with the visual aspect of Weinberg & Schmeling (1992) is not as good as for cases I and II. The intermediate fluid coalesces at the bottom of the rising diapir instead of being carried upwards as in Weinberg & Schmeling (1992). However, the wavelength of interface 1–2 is similar. Using a more refined grid or an arithmetic interpolation for the viscosity at the interface did not improve the results and cannot explain the observed discrepancy. In fact, we obtain similar results when we use a viscosity ratio γ = 100 and set random perturbation only on the lower interface. Explanations for this discrepancy can be due either to i) erroneous setup details provided in Weinberg & Schmeling (1992) or ii) instabilities tending to develop faster in JADIM and OpenFOAM than in Weinberg & Schmeling (1992).

3.3.3 Comparison with RT linear stability theory

In the following, we focus on the initial development of the instabilities at interfaces 2–3 and 1–2. Their growth rates |$\sigma _{23} = \log (V_{2-3}/\mathcal {U})/t$|⁠, |$\sigma _{12} = \log (V_{1-2}/\mathcal {U})/t$| and their wavelengths λ₂₃ and λ₁₂ are compared with those given by the linear stability theory developed by Ramberg (1981a) (cf. section 3.1 and Appendix A1) for three-layer systems. The dominant wavelengths λ₂₃ and λ₁₂ for each interface 2–3 and 1–2 both correspond to the maximum growth rate.

In the numerical models, the wavelength is measured as the distance separating the axes of symmetry of the diapirs. In case I for instance, one observes at time |$t/\mathcal {T}\approx 10$| that the selected wavelength for JADIM is λ/h₂ ≈ 2.1, while for OpenFOAM, one finds λ/h₂ ≈ 2.7 (Fig. 4a). The growth rate is obtained from picking at each time step the velocity of the maximum vertical location of each interface 1–2 and 2–3. The evolution of interface 1–2 is displayed in Figs 5(a), (c) (e) for cases I, II and III, respectively. The evolution of interface 2–3 is similar to that of interface 1–2 and it is not shown.

Figure 5.

Three-phase Rayleigh–Taylor system by Weinberg & Schmeling (1992). Left-hand column: time evolution of the velocity of the maximum height of interface 1-2, simulated in Fig. 4 (harmonic viscosity averaging is used here). Right-hand column: maximum growth rate obtained in the simulations as a function of the wavelength, for interfaces 1-2 and 2-3. (a-b) Case I, (c-d) Case II, (e-f) Case III. The red curves correspond to the linear stability theory (Ramberg 1981a).

The lin-log representation allows to straightforwardly observe any exponential growth (curve with a constant slope) and measure the growth rate. In case I for instance, one observes a first instability at |$5 \lesssim t/\mathcal {T} \lesssim 20$| (Fig. 5a). A second exponential growth is seen for |$50 \lesssim t/\mathcal {T} \lesssim 200$|⁠. This is the signature of the double overturn: the first slope describes when layer 2 grows through layer 3, the second slope when layer 2 grows through layer 1.

For case II, JADIM and OpenFOAM present very similar growth rates and wavelengths until t/τ = 100. Then, ASPECT (orange curve) exhibits a higher growth rate than OpenFOAM (green curve) and JADIM (blue curve). For comparison, we plot in the left column of Fig. 5 the prediction of the linear stability theory (red and pink lines for interfaces 1–2 and 2–3, respectively). The slopes of the curves are important whereas the initial location at the origin is not. If we plot parallel lines to the red curves, all codes are in good agreement. However, we find that ASPECT presents some oscillations at t/τ < 150, that is at the onset of the destabilization. Therefore, for this code, we consider the slope of interface 2-3 between time 70 < t/τ < 200, and the slope of interface 1–2 at 150 < t/τ < 370, to evaluate the growth rate displayed in Fig. 5(d).

In order to provide a more quantitative comparison, Figs 5b, d, f displays the dispersion relation (growth rate versus wavelength) obtained from the linear stability theory of Ramberg (1981a) (red curves). For cases II and III, one curve displays the two main wavelengths corresponding to each interface 1-2 and 2-3. For case I, two curves are needed since only interface 2-3 deforms and the instabilities occur at two different times (one time for layer 3 to cross layer 2 and one time for layer 3 to cross layer 1). The growth rates and wavelengths in our simulations are also plotted (symbols). It is worth noting that measuring such growth rates from the simulation is somewhat difficult as its specific value depends on the time span chosen for computing the slope (see Fig. 5, left column) and we have chosen to use the maximum growth rate of each interface. In any case, this at least partly explains why the values obtained in the simulations can be slightly larger (or smaller) than the maximum growth rate predicted by the linear stability theory.

In summary, the agreement between the wavelength and growth rate obtained in the simulations using our VOF codes and those predicted by the linear stability analysis of Ramberg (1981a) is rather good since it is always within 15 − 20 per cent. ASPECT does not produce such a good fit (see the growth rate of interface 1-2 in Fig. 5d, triangle). The discrepancies can be mainly explained by the difficulty to accurately measure the growth rate since it varies over time. Once again, even if the visual aspect of case III is different from that of Weinberg & Schmeling (1992) (Fig. 4c), the growth rate and wavelength given by JADIM and OpenFOAM remain close to the linear stability prediction.

For completeness, we compare our three-layer system results with alternative theories mentioned in the introduction (Ramberg 1981a; Burg et al. 2004), and therefore separate our system into a top layer and a bottom layer. There are then two ways to proceed: either group the two lower layers together (top layer = layer 1 and bottom layer = layers 2 and 3) or, since the viscosity ratio between layers 1 and 2 is high, layer 1 is set to act as a wall, layer 2 is considered as the top layer and layer 3 as the bottom layer. Hence we compare the growth rate and the wavelength obtained with our codes (JADIM, OpenFOAM, ASPECT) for case II, with:

• Burg et al. (2004)’s predictions for a two-layer system (eq. A2, referred to as BKP),

• Ramberg (1981a)’s theory for a two and a three-layer system (eq. A1, referred to as R).

Results are presented in Table 2 a,b for interfaces 1-2 and 2-3, respectively. Note that for interface 1-2, we consider layer 1 as the top layer, and layers 2 and 3 as the bottom layer; for interface 2-3, we consider layer 2 as the top layer and layer 3 as the bottom layer. The wavelength of interface 2-3 for a three-layer system is correctly predicted by the two-layer system theory of Ramberg (1981a), but the growth rate and wavelength of interface 1-2 are very different. BKP’s predictions of both the wavelength and growth rate are very different from the numerical results and do not seem appropriate to describe the present three-layer system. It turns out that the multi-layer theory of Ramberg (1981a) is the one found able to predict both the wavelength and the growth rate for a three-layer system. Indeed, results given by our codes and those from Ramberg (1981a)’s theory for a three-layer system are in good agreement for both interfaces.

Table 2.

Theoretical and modelled maximum growth rate κ₁₂ and corresponding wavelength λ₁₂ of interface 1-2 (a) and 2-3 (b) Rayleigh–Taylor instabilities case II; BKP: Burg et al. (2004); R: Ramberg (1981a).

(a) Interface 1–2
	JADIM		OpenFOAM		ASPECT		2 layers theory		3 layers theory
	2-D	3-D	2-D	3-D	2-D	3-D	BKP	R	R
λ12/h2	13.2	> 9	14.4	> 9	11.32	> 9	7.2	42.0	15.2
κ12	0.017	0.0165	0.017	0.015	0.018	0.02	0.98	13.7	0.017
(b) Interface 2-3
	JADIM		OpenFOAM		ASPECT		2 layers theory		3 layers theory
	2-D	3-D	2-D	3-D	2-D	3-D	BKP	R	R
λ23/h2	2.4	2.4	2.6	2.4	2.0	2.24	0.8	2.6	2.6
κ23	0.015	0.015	0.016	0.0158	0.0148	0.0164	0.1	0.153	0.0154

3.3.4 Numerical sensitivity of the results - RT case II

Here we carry out various sensitivity tests for case II with JADIM, OpenFOAM and ASPECT. We visually compare structures as well as the evolution of velocity over time. We plot the mass error, and assess the influence of i) the mesh resolution, ii) the initial perturbation of the interface and iii) the viscosity averaging at the interface (harmonic or arithmetic).

• Mass Error - The time evolution of the relative mass error |$\displaystyle {\Delta M^i = (M^i-M_0^i)/(M_0^i)}$| with Mⁱ the mass of phase i and |$M_0^i$| the initial mass of phase i, is presented in Fig. G1. Both JADIM and OpenFOAM display values lower than 10⁻⁷.

• Mesh resolution - Grid convergence is illustrated in Table E1 of Appendix  E. Grid independent results are obtained with a relatively coarse grid for JADIM and OpenFOAM (90 × 37).

• Initial Perturbation of the interfaces - We plot and discuss with figures in Appendix  E the time evolution of the velocity of the maximum vertical location of interfaces 1-2 and 2-3, with or without perturbing the two interfaces with a random perturbation in JADIM. In summary, simulations with an initial perturbation are in better agreement with the theory (Figs 5 c–d). In contrast for cases I (Fig. E1a) and III (Fig. E1c)) the shape and dynamics of the interfaces is found to evolve rather independently from initial perturbations at the interfaces.

• Influence of the type of viscosity averaging - Fig. E3 of Appendix  E compares the time evolution of the velocity of the maximum vertical location of interface 1-2 for an arithmetic and a harmonic viscosity averaging at the interfaces. The evolution of interface shapes are roughly similar, yet, harmonic averaging favours a behaviour controlled by the lowest viscosity at an interface, rendering it easier for layer 2 to invade layer 1 and the corresponding diapir becoming bigger, in comparison with results produced with a arithmetic viscosity averaging. However, we note that arithmetic averaging provides results more similar to Weinberg & Schmeling (1992) at t/τ = 497.

3.3.5 Comparison between 2-D and 3-D RT models

In this section, we perform 3-D simulations extending the 2-D setup of Weinberg & Schmeling (1992) for cases I, II and III. The first objective is to show the capability of our VOF methods to tackle 3-D Rayleigh–Taylor problems; the second objective is to assess the sensitivity of the theoretical predictions to possible 3-D effects. The initial and boundary conditions are similar to the 2-D cases. The size of the domain however, has been shortened along the horizontal directions and L = 1.12H and is discretized using a resolution of 112 × 112 × 100.

The temporal evolution of the 3-D system of cases I, II and III are displayed for OpenFOAM in Fig. 6. Case II is displayed in Fig. 7 for JADIM and ASPECT. In all cases, layer 3 starts to destabilize by forming small domes which later evolve as individual diapirs, as previously shown by Biot (1966); Berner et al. (1972); Talbot et al. (1991); Kaus (2004); Fernandez & Kaus (2015).

Figure 6.

3-D Rayleigh–Taylor simulations of cases (a) I, (c) II and (e) III corresponding to the 2-D analogues presented in Fig. 5 (Weinberg (1992)) using OpenFOAM with harmonic viscosity averaging. Here, the horizontal size of the domain is 1.12H and the grid size is 112 × 112 × 100. Temporal offsets are chosen, in 3-D, for case I: |$\text{Off}_{OP}^1 = 0$|⁠, case II: |$\text{Off}_{OP}^2 = 300$|⁠, case III: |$\text{Off}_{OP}^3 = 216$|⁠. The offsets in 2-D are, for case I: |$\text{Off}_{OP}^1 = 5$|⁠, case II: |$\text{Off}_{OP}^2 = 276$|⁠, case III: |$\text{Off}_{OP}^3 = 300$|⁠.

Figure 7.

3-D Rayleigh–Taylor simulations of case II corresponding to the 2-D analogues presented in Fig. 5 (Weinberg (1992)) using JADIM (a-b) and ASPECT (c-d) with harmonic viscosity averaging. Here, the horizontal size of the domain is 1.12H and the grid size is 112 × 112 × 100. Temporal offsets are |$\text{Off}_{JA}^2 = 282$| in 3-D and |$\text{Off}_{JA}^2 = 300$| in 2-D.

In case I (Fig. 6a), since layer 2 is denser than the other layers, layer 1 acts as a boundary and the little diapirs merge below interface 1-2 before rising through layer 1. In case II (Figs 6 d and 7b, d), the little diapirs grow without merging and layer 2 entrains the diapirs on a side of the box (as in 2-D). Then, both layers 2 and 3 rise through layer 1. In case III (Fig. 6e), the dynamics is roughly similar to that of case II with a larger delay due to the larger viscosity ratio between layers 1 and 2 (200). The little diapirs merge and rise through layer 1 at the same time.

For case II, JADIM, OpenFOAM and ASPECT present similar patterns of deformation, and the general dynamics is quite similar in 2-D and 3-D. In addition, the wavelength of the patterns are of the same order of magnitude, as indicated in Tables 2 and 2. Biot (1966) extended to 3-D his 2-D linear stability solutions for a multi-layer system. He showed that the distance between peaks and crests in 3-D is almost that of the wavelength of the 2-D solution (λ_2-D) as observed here. The spacing of diapirs both in 2-D and 3-D follows the predictions of the linear stability theory (if they grow from an initial interface perturbed with a random perturbation).

The time evolution of the interfaces velocities in 3-D are compared with the 2-D results in the left column of Figs 7 and 6. Interestingly, this figure shows that the first instability grows at a similar rate in 2-D and in 3-D (slopes are similar). However, the 2-D and 3-D dynamics are not strictly identical: in 3-D the rising velocity of the second instability may be different from the 2-D one (when fluid 3 has already intruded fluid 2, cf. also second exponential in cases I and III), as shown by the offset between curves. The observed relative difference in velocity can be as large as 15 per cent. Moreover, if JADIM and OpenFOAM yield continuous curves, this is not the case for ASPECT, where, at time 100 < t/τ < 200, curves are shifted respect to their continuation.

Fig. 6 b displays the pattern of layer 3 inside layer 2 for case I (note that we do not consider the pattern of layer 3 into layer 1 since the wavelength of the corresponding instability is of the order of the box size): a polygonal pattern is observed. Biot (1966) studied 3-D patterns of a two-layer Rayleigh–Taylor instability and found that for triangular or hexagonal patterns, the theoretical distance h separating two neighbouring peaks in one horizontal direction relative to that in the orthogonal direction l follows the relation |$\displaystyle {h/l= 1.155}$|⁠. Here, we measure the average distances |$\displaystyle {h^m/h_2\approx 2.48}$| and |$\displaystyle {l^m/h_2\approx 2.16}$|⁠, thus giving |$\displaystyle {h^m/l^m \approx 1.148}$| which is close to the value of Biot (1966). In addition, Biot (1966) states that the characteristic distance d between peaks is determined to be d = 1.155λ_2-D. Here, by computing the mean distance between peaks, we recover that value.

As found in the previous section, JADIM and OpenFOAM produce equivalent results compared to ASPECT. Here again, the mass error (in per cent) for JADIM and OpenFOAM remains below 10⁻⁷ in 3-D whereas it is only below 10⁻³ for ASPECT (Fig. G1).

4 MULTILAYER RAYLEIGH–BÉNARD INSTABILITIES

In this section, we consider two configurations involving Rayleigh–Bénard instabilities: the case of a single layer of fluid heated at its base with a punctual source, to be compared to the work of Vatteville et al. (2009), and the case of an initially stable stratified two-layer system heated from below, to be compared to the work of Le Bars & Davaille (2004).

In order to solve such configurations with our VOFs methods, the full set of eqs (2)–(5) is solved and both the density and viscosity are now prescribed temperature-dependent. In particular, the density follows a Boussinesq law |$\tilde{\rho }= \rho _0 \cdot \mathcal {F} (T)$| with |$\mathcal {F}(T) =1-\alpha (T-T_{ref})$| where T_ref is a reference temperature and α is the thermal expansion coefficient. The viscosity follows a law |$\tilde{\mu }= \mu _0 \cdot \mathcal {G} (T)$|⁠, which will be detailed for each configurations.

4.1 One-layer Rayleigh–Bénard system

We aim here at reproducing the well-controlled laboratory and numerical experiments performed by Vatteville et al. (2009). Initially, a fluid of density ρ₁,viscosity μ₁ and temperature T_in is placed in a cylindrical tank of radius L and height H (see Fig. 8a). A local heat source is placed in a small region (r_h) at the centre of the cylinder’s base, with a temperature that increases with time from T_in to T_H following the law T(t)/T_R = 1.0 + 0.09(1 − exp ( − t/26.3)) (T in Kelvin). This problem is controlled by two dimensionless parameters, for example the Prandtl number, a function of ν = μ/ρ the kinematic viscosity and D_th the heat diffusivity, and the Rayleigh number, where ΔT = T_H − T_in:

$$\begin{eqnarray*} { Pr=\frac{\nu }{D_{th}}}, \quad { Ra = \frac{\alpha g \Delta T H^3 }{ D_{th} \nu }}. \end{eqnarray*}$$

(12)

Figure 8.

Rayleigh–Bénard problem (Vatteville et al. 2009): comparisons of the velocity field between laboratory (Vatt-Exp) and numerical experiments (Vatt-num), OpenFOAM and JADIM. Velocity is scaled with V_Stokes = αgΔTH²/ν_max, time is scaled with τ = ν_max/(αgΔTH), and distance is scaled with the domain height H. The grid size in JADIM and OpenFOAM is 170 × 322.

The problem is assumed to be axisymmetric around a central vertical axis, therefore computations are performed using a 2-D grid of size 170 × 322 along the radial and vertical directions, respectively. Boundary conditions are given in Fig. 8(a). The temperature-dependent viscosity |$\tilde{\mu }= \mu _0 \cdot \mathcal {G} (T)$| is provided by Vatteville et al. (2009), who evaluated an empirical law |$\mathcal {G} (T) = 1.9 \exp (-7.11+1892.0/T)$|⁠, T in Kelvin. With their choice of temperature range [T_in, T_H], the largest viscosity contrast reaches μ_max/μ_min ≤ 2. Grid convergence is illustrated for JADIM in Fig. F1.

In the experiment, and after some time during which heat diffuses in the thermal boundary layer around the heat source, a plume grows until it reaches a diameter of the order of that of the heat source, and rises. Fig. 8(b) displays the temporal evolution of the maximum velocity along the axis of symmetry, denoted ’conduit’, in Vatteville et al. (2009)’s experiments and in the JADIM and OpenFOAM tests. Velocity isocontours are plotted over the whole model domain at a specific time (Fig. 8c), and then along the axis of symmetry at three different times (Fig. 8d).

The ‘shape’ of the plume (defined with the iso-values of velocity in Fig. 8(c) obtained with both our codes compares well with that of Vatteville et al. (2009)’s experiment. The maximum velocity along the plume conduit as a function of time is presented in Fig. 8(b). The figure shows, for JADIM and OpenFOAM, a maximum velocity higher than in the lab experiments but lower than in the numerical simulation of Vatteville et al. (2009). At later times, good agreement is found for all numerical approaches. Their value however is larger than the one measured in the experiment (by about 6 per cent). This discrepancy was attributed by Vatteville et al. (2009) to the laboratory measurements which made use of a local spatial averaging procedure which tends to moderately underpredict the velocity.

4.2 Two-layer Rayleigh–Bénard system

We consider now the more complex problem of a two-layer Rayleigh–Bénard configuration. In this section, we set up numerical cases comparable to the laboratory experiments of Le Bars & Davaille (2004). Although this work was oriented to model the Earth’s mantle, it may very well constitute a good basis to explore the conditions for crustal scale convection, as will be discussed later.

4.2.1 Case description

A square tank of horizontal width L and height H is filled with two fluid layers, the uppermost layer of density ρ₁, viscosity μ₁ and height h₁ being less dense and less viscous than the bottom layer of density ρ₂, viscosity μ₂ and height h₂ (Fig. 9a). The system is initially at a constant temperature T_in, and is simultaneously heated from below at temperature T_H and cooled from above at temperature T_c. Note that here the system is initially stable. Assuming that the thermal expansion coefficient α and the heat diffusivity D_th are the same in both fluids, the present problem is controlled by five dimensionless parameters, namely the Prandtl and the Rayleigh numbers (eq. 12), the thickness and viscosity ratios:

$$\begin{eqnarray*} {a=h_2/H, \gamma =\mu _1/\mu _2}, \end{eqnarray*}$$

(13)

and the buoyancy number:

$$\begin{eqnarray*} {B = \frac{\rho _2 - \rho _1}{\alpha \Delta T(\rho _1+\rho _2)/2 }}, \end{eqnarray*}$$

(14)

where ΔT = T_H − T_C. The parameter B expresses the ratio of the stabilizing density stratification to the destabilizing thermal density anomaly. Here, time is scaled by τ_T = μ₂/(αgΔTH) and the origin of time t = 0 has been set as the time when the interface height y reaches y/h₂ = 1.08. Le Bars & Davaille (2004) performed various experiments which enabled them to draw a phase-map of the various dynamics which may occur, depending on a, B and, to a lesser extent γ and Ra (Fig. 9b). In particular, one can see that B controls the ability of the interface to deform: if B > 0.3, the interface remains mostly flat, and convection occurs in both layers separately. If B ≤ 0.3, the interface is fully destabilized and convection occurs in the whole domain. Their parameters vary between: 0.08 < a < 0.95, 7 · 10⁻⁴ < γ < 1, 10³ < Ra < 10⁸, 0.1 < B < 4, Pr > 100. Below, we model four cases with fixed a = 0.3, γ ≈ 10⁻³, Pr ≈ 10⁶ and Ra ≈ 10⁵ and vary the buoyancy number in the range 0.2 ≤ B ≤ 2.

Figure 9.

The two-layer Rayleigh–Bénard problem by Le Bars & Davaille (2004): (a) Initial setup, geometry and boundary conditions (L/H = 2. μ₀ = 0.5 − 5 Pa.s and |$\mu = \mu _0 \cdot \mathcal {G} (T)$| with |$\mathcal {G} (T) = 2.2 \cdot e^{-0.038T}$|⁠, T in Celsius). (b) Regime diagram of the convection regimes according to the lab experiments in the parameter space (B; a), modified from Fig. 2 of Le Bars & Davaille (2004). The coloured stars indicate the location of the 3-D simulations performed in this work. (c) Distribution of the instantaneous temperature field in a vertical plane at time |$t/\mathcal {T} = 1800$| computed with JADIM (γ = μ₁/μ₂ ≈ 10⁻², Ra = 1.6 × 10⁵, Pr ≈ 3.5 × 10⁶). The grid size in JADIM and OpenFOAM is 90 × 90 × 44.

4.2.2 Qualitative results

3-D simulations are performed with the geometry, initial and boundary conditions displayed in Fig. 9(a) and described above. The computational domain is discretized using a 90 × 90 × 44 grid-size. Following Le Bars & Davaille (2004) we use an exponential temperature-dependent viscosity such that μ = μ_o · exp ( − 0.038T).

We set four simulations corresponding to different locations in the (B, a) parameter space (see star symbols in Fig. 9b): B = 0.2, 0.3, 0.6 and 2. Snapshots of the typical modelled structures are illustrated Fig. 9(c). Our simulations produce deformation modes in agreement with those identified by Le Bars & Davaille (2004), in particular:

• Case I, B = 0.2 (equivalent to ρ₂ < ρ₁, Fig. 10a): the density of the fluid from below decreases enough to invade the fluid from above as a diapir. At the end of the experiment, a whole single-layer convection mode develops.

• Case II, B = 0.3 (equivalent to ρ₂ ∼ ρ₁, Fig. 10c): the system stands at the transition between the interface remaining flat or deforming like a rising diapir. Here, the fluid from below intrudes the fluid from above, in a diapiric manner.

• Case III, B = 0.6 (equivalent to ρ₂ > ρ₁, Fig. 10b): both layers remain stagnant but the interface progressively deforms as convective structures develop in each layers.

• Case IV, B = 2 (equivalent to ρ₁ ≪ ρ₂, Fig. 10d): convection appears separately in each layer and the interface remains flat.

Figure 10.

Rayleigh–Bénard two-phase system (Le Bars & Davaille (2004)): Iso-contour C = 0.5 for various B: 0.2(a), 0.3(b), 0.6(c) and 2.0(d). Time is scaled by τ = ν₂(T_H)/(αgΔTH).

Fig. 10 displays the 3-D flow dynamics for each for the four tested B cases. JADIM and OpenFOAM behave similarly: both cases I (Fig. 10a) and II (Fig. 10b) with B ≤ 0.3, show deformation modes with diapirism followed by whole-layer convection. While JADIM diapirs rise directly to the top (Fig. 9c), OpenFOAM displays diapirs that collapse on themselves before rising all the way to the top (see Figs 10a and b, time step 8200).

ASPECT is also tested for case II (Fig. 10b). Thermal plumes form and fall from the top, deforming the interface before diapirs can grow from layer 2. Then, as with JADIM, diapirs grow and reach directly the top of the box. A large random perturbation (2/3 of a cell size) has to be inserted in order to obtain a similar development of diapirs to JADIM or OpenFOAM. Otherwise (with a smaller random perturbation), a single diapir grows from the centre of the box, later surrounded by five emerging diapirs. This different behaviour of ASPECT is also discussed in comparison with analytical predictions in the next section.

Cases III and IV with B > 0.3 display an interface that only slightly deforms over time as expected from Le Bars & Davaille (2004, Figs 10c, d and 9c).

4.2.3 Comparison with analytical predictions

In cases I and II, deformation of the fluid interface is significant, as layer 2 progressively becomes less dense than layer 1. At this point, the system becomes similar to a Rayleigh–Taylor system. Thus, we may compare the modelled growth rate of the interface to linear stability analytical solutions (LST, see Section 3.1), would not there be the difficulty rising from temperature-dependent densities and viscosities. Since we haven’t found analytical solutions for exponentially temperature-dependent viscosity (closest predictions might be Popov et al. (2014) for exponentially depth-dependent viscosity), we choose to first compare numerical growth rates in between codes, and then to identify the associated range of temperature with which the latter can be delimited by Ramberg (1981a)’s analytical growth rate.

We thus proceed for case II (B=0.3), and display in Fig. 11 the modelled and theoretical maximum heights of interface 1–2 with time. The evolution of this interface can be divided in two stages, before and after t/τ ≃ 2500 (see Figure caption for the choice of time origin). Before that time, the system does not significantly destabilize, and the interface 1–2 remains relatively flat with JADIM an OpenFOAM. In contrast with ASPECT, this interface deforms more, and thermal plumes are seen to drip from the top before any instability actually manages to grow there. After that time, diapirs rise, first in JADIM (from t/τ ≃ 2250), then in OpenFOAM(from t/τ ≃ 2500), and finally in ASPECT (from t/τ ≃ 2800). From then on the interface growth rates are close to each other by ∼10 per cent, namely κ_JA = 0.10 for JADIM, κ_OP = 0.12 for OpenFOAM and κ_AS = 0.16 for ASPECT.

Figure 11.

Two-layer Rayleigh–Bénard system from Le Bars & Davaille (2004), B = 0.3, case II. Temporal evolution of the highest position of the layers’ interface for JADIM, OpenFOAM and ASPECT. Δy/Δy_max = (y(t) − y₀)/(H − y₀), y₀ being the interface’s initial height and H the box’s total height. The origin of time t/τ = 0 has been chosen so that the interface reaches y/h₂ = 1.2 at the same time for all codes. Dotted lines are the tangents to the modelled curves. Red lines referred as LST provide the theoretical growth rate slopes (Ramberg 1981b) for |$T_{R}=93\% \Delta T$| and T_R = ΔT with ΔT = T_H − T₁. The two orange curves for ASPECT correspond to a weak random perturbation ‘wwn’ (6 per cent of the cell size) and a strong random perturbation ‘swn’ (66 per cent of the cell size), showing both high initial amplitudes prior to the onset of destabilization.

With ASPECT, we tried to diminish the amplitudes related to the initial destabilization stage (t/τ ≤ 2800), by reducing the amplitude of the initial random perturbation: a single diapir is seen to form at the centre of the model domain prior to others, and then the transition to the main growth rate is shifted by t/τ = 1000 (compare orange curves in Fig. 11). Since on the other hand, mass error in ASPECT is found to be about 10⁻² (Fig. G1c), we believe that diffusion may be causing this less ’stable’ initial growth compared to the two other VOF codes. Further tuning of specific numerical settings in ASPECT might improve its behaviour, yet we did not pursue this matter further as this was not the aim of our study.

Then we plot in Fig. 11 the LST growth rates that delimit the slopes of our modelled growth rates. Since this LST growth rate depends on ρ₁, μ₁, ρ₂, μ₂ (eq. 10 and Appendix  C), which all vary with temperature, we need to choose an equivalent temperature for layer 1 and an equivalent temperature T_R for layer 2. For layer 1, we take |$T_1 = 34 \% (T_H-T_C)$| (which is its average temperature before destabilization) to deduce ρ₁ and μ₁. For layer 2, we seek the temperature T_R (which determines ρ₂, μ₂) with which we can trace the lower and upper bounding slopes of our modelled growth rates: we find |$\kappa _{T_R}= 0.09$| and |$\kappa _{T_R}= 0.18$|⁠, associated to reference temperatures |$T_{R}=93\,{\rm{per\,cent}} \cdot (T_H-T_1)$| and |$T_{R}=100\,{\rm{per\,cent}} \cdot (T_H-T_1)$|⁠, respectively. This allows us to conclude that the numerical codes reproduce the development of RB instabilities with a precision of 7 per cent of the ‘equivalent stratified’ theoretical growth rate. We cannot state more precisely which code best matches a true solution, but at least JADIM, OpenFOAM and ASPECT present an overall consistent evolution of their interface, providing confidence in their behaviour.

Concerning the dominant wavelength λ, the linear stability theory of Ramberg (1981a) predicts that λ/h₂ = 3 for a Rayleigh–Taylor system. We obtain for OpenFOAM and ASPECT λ/h₂ = 1.5, and for JADIM λ/h₂ = 1.2. However, theoretical values for Rayleigh–Taylor instabilities do not correspond to those for a two-layer Rayleigh–Bénard system, as shown by Le Bars (2003). In fact Le Bars (2003) determined experimentally that λ/h₂ = 9.1 × Ra^−0.14 which in our case leads to λ/h₂ = 1.7. OpenFOAM and ASPECT provide a value that differs by 12 per cent to this experimental law, while JADIM differs by 30 per cent.

5 SYNTHESIS OF CODES PERFORMANCES

We have produced numerical models with two VOF codes, JADIM and OpenFOAM, for two- and three-layer systems with and without thermal effects, and have shown a good agreement with the previous studies of van Keken et al. (1997), Weinberg & Schmeling (1992), Ramberg (1981a), Vatteville et al. (2009) and Le Bars & Davaille (2004). The two codes present some differences in their implementation, in particular the treatment of the transport equation of volume fraction (eq. 2). For instance, in JADIM, the ‘numerical thickness’ of the interface is larger than that in OpenFOAM (2–3 grid cells versus 1 grid cell) and this may explain the slight differences observed in the numerical results between both codes (e.g. Figs 2 and 10). If we compare the results obtained with both VOF methods and those obtained with a field method like ASPECT, the field method tends to be more diffusive (>3 grid cells for the interface).

We report below the technical performances of JADIM, OpenFOAM and ASPECT: their weak and strong scalabilities as well as the computational time required for some of the experiments described above.

5.1 Scaling

5.1.1 Strong scaling

We assess the strong scalability of JADIM, OpenFOAM and ASPECT in the case of the simulation displayed in Fig. 10(b), case II of the two-layer Rayleigh–Bénard problem of Le Bars & Davaille (2004). More precisely, we fix the size of the computational domain to 120 × 120 × 60 gridpoints (i.e. about 864 000 gridpoints in total), and measure the computational time with respect to the number of processors, denoted N_p. Fig. H1 a displays the computational time required to do ten iterations as a function of N_p. A perfect scaling would lead to a line log–log of slope 1/N_p (Fig. H1a, dotted line). The computational time is scaled by the time required using 4 processors (two processors on two different nodes) since the available memory on one node (≈192 Go) was not sufficient to support a run with ASPECT.

For N_p ≥ 8, the speedup decreases for JADIM, OpenFOAM and ASPECT. This indicates that when the equivalent problem size per processor is about 50³ or less, the time devoted to communication between processors is of the same order of magnitude or larger than that devoted to solve the equations. Using a larger domain size, involving 10⁶ or 10⁷ gridpoints, improves the performance presented in Fig. H1(a). Furthermore, OpenFOAM and ASPECT speed up twice faster than JADIM at large N_p, and as N_p ≥ 32 (i.e. 30³ gridpoints per processor), with an advantage for OpenFOAM.

5.1.2 Weak scaling

We assess the weak scalability of JADIM, OpenFOAM and ASPECT using again the case displayed in Fig. 10(b) [case II of the two-layer Rayleigh–Bénard problem of Le Bars & Davaille (2004)]. We fix the grid size of the computational domain per processor to 30 × 30 × 30, and vary the computational domain size with the number of processors. We measure the solution time with respect to the number of grid cells. Fig. H1(b) displays the computational time required to achieve ten iterations (scaled with the computational time on one processor). A perfect scaling would lead to a horizontal line. For both codes, the parallel efficiency decreases as we increase the number of processors due to communication between processors. For 108 000 grid cells (i.e. four times more cells than in the case of Section 4), 50 per cent of the time is used for communication in JADIM, 20 per cent for OpenFOAM and 10 per cent for ASPECT. For 32 processors, the speedup loss is about 40 per cent for ASPECT and 90 per cent for JADIM.

5.2 Computational time in 3-D

We compared computational times between codes for different test runs and found that OpenFOAM is in general faster than JADIM (see e.g. Fig. H1) and ASPECT. For instance, the 3-D computations of Weinberg & Schmeling (1992)’s configuration presented in Figs 7 and 6 took, for the same physical time and the same time step, on Intel(R) IVYBRIDGE 2.8 GHz processors:

• 6 days on 8 cores with OpenFOAM versus 14 days on 27 cores with JADIM (case I)

• 8 days on 8 cores with OpenFOAM versus 24 days on 27 cores with JADIM versus 8 d on 18 cores with ASPECT (case II)

• 8 days on 8 cores with OpenFOAM versus 18 days on 27 cores with JADIM (case III)

In this case, ASPECT needed 18 cores to achieve the same physical time as OpenFOAM. JADIM needed much more time and cores than OpenFOAM and ASPECT. Note that even if JADIM and OpenFOAM use the same numerical scheme to compute the viscous term, that of JADIM is not parallel, in contrast to OpenFOAM. In addition, we find that JADIM’s pre-conditioner used to solve pressure (conjugate gradient with block Jacobi method) implies that with an increasing number of cores the number of iterations for convergence increases. This is not the case for OpenFOAM which uses a stabilized preconditioned biconjugate gradient with a Diagonal-Incomplete Cholesky pre-conditioner. These two points may explain the poor performances of JADIM.

Additionally, the 3-D computations of Le Bars & Davaille (2004)’s case II (B=0.3, Fig. 10) took:

• 37 days on 4 cores with JADIM

• 13 hours on 4 cores with OpenFOAM

• 15 days on 32 cores with ASPECT

Computational times here are only informative since (i) all processors were different and (ii) we realize that some numerical settings might have been better optimized. For instance we note that the time step of JADIM computations was set limited 10 times smaller than that of OpenFOAM. Together with the result from Fig. G1 showing that JADIM is about two times less efficient than OpenFOAM, the different computational times obtained here in 3-D are justified.

Furthermore in 3-D here, OpenFOAM is found to run faster than ASPECT, despite we obtained a better parallel efficiency for ASPECT in the previous section. Together with the peculiar results of this ASPECT case (large initial instabilities preceding diapir rise, Section 4.2.3), we conclude that some memory allocation and computational options of ASPECT aught to be further investigated.

5.3 Synthesis

In summary, ASPECT runs faster for 2-D simulations than JADIM and OpenFOAM. But, JADIM and OpenFOAM better conserve mass than ASPECT.

OpenFOAM seems the best adapted code to study crustal polydiapirism since (i) it is faster than ASPECT and JADIM and ASPECT in 3-D and (ii) it conserves mass better than ASPECT. In the following section, we thus use OpenFOAM to model Naxos’s observations.

6 APPLICATION TO THE DEVELOPMENT OF METAMORPHIC DOMES IN NAXOS

Vanderhaeghe et al. (2018) interpreted the domes of Naxos Island, Greece, as the exhumation of imbricated or adjacent polydiapirs in a larger rising dome. According to Vanderhaeghe et al. (2018), these domes would have formed after crustal thickening dated at 55 Myr, and comprised a first episode of crustal scale convection from 24 to 16 Myr, followed by a second episode of polydiapirism from 16 to 13 Myr, possibly associated with thinning of the orogenic crust. An estimate of the characteristic size and growth rate of these diapiric instabilities was proposed by Vanderhaeghe et al. (2018) using the critical Rayleigh number threshold; a ‘convectable’ crustal thickness H between 10 and 30 km requires a viscosity range of 10¹⁶−10¹⁸ Pa.s and a density contrast of 50−200 kg m⁻³. Vanderhaeghe et al. (2018) estimated that the large dome covers an elliptic area of dimension L_Na = 5 × 12 km, that the subdomes have a size d_Na = 2−3 km, and that the velocity of convection was v_Na = 1−5 cm yr⁻¹. A revolution period of about 2 Myr of convective cycles was estimated based on zircon geochronology (Vanderhaeghe et al. 2018). More precisely, two superimposed destabilization processes may have been at play:

• A convection episode, analog to a destabilizing single layer Rayleigh–Bénard system with thermally dependent density and viscosity,

• A polydiapirism episode, analog to a case-II three-layer system defined by Weinberg & Schmeling (1992), Fig. 6(c).

The increase in temperature responsible for these gravitational instabilities may have been caused by an increasing mantle heat flux and/or the thickening of a radiogenic crust (Thompson & England 1984; England & Thompson 1986; Ueda et al. 2012). The models proposed here present both mechanisms of convection and polydiapirism independently. Our aim is to better constrain the propositions of Vanderhaeghe et al. (2018) for the case of Naxos, with the support of 3-D models using OpenFOAM. Although our results already provide some insight on the dynamics of Naxos’s crust, another study will be required in order to explore and demonstrate how the transition or the combination of both these episodes would have occurred.

6.1 Naxos Episode 1, Rayleigh–Bénard convection

For the convection episode, we build a single-layer Rayleigh–Bénard system representing a pre-thickened orogenic crust. We assume that the first 10 km of cold upper crust are not affected by deeper mechanisms of deformation, therefore the model domain starts 10 km below the Earth’s surface. The modelled crustal thickness is either 35 or 50 km, with an initial linear geotherm spanning T_C = 300 °C at the upper boundary (–10 km depth), and T_H = 600 °C at the bottom boundary, down to depth –45 or –60 km. Kinematic boundary conditions are prescribed: no slip at the lower boundary and free slip at the lateral and upper sides.

The densities are defined according to the Boussinesq approximation (eq. 7), with a reference density |$\rho _1^0$| given at T_C = 300 °C. The viscosities are defined bounded by maximum and minimum values related to minimum and maximum temperatures T_C and TH, according to an Arrhenius law |$\mu _i = \displaystyle {B \cdot e^{\frac{A}{T}}}$|⁠, with |$\displaystyle {A = \frac{\ln (\mu _{\rm{max}}/\mu _{\rm{min}})}{1/T_{\rm C}-1/T_{\rm H}}}$| and |$\displaystyle {B = \frac{\mu _{\rm{min}}}{e^{A/T_{\rm H}}}}$|⁠. A and B are ajusted so that viscosities vary in the range 10¹⁷ Pa.s (at 900 or 1000 °C)–10²¹ Pa.s (at 300 °C). Depth-dependent profiles of the initial and final temperature, the density and the viscosity are presented in Figs 13(a)–(c).

At time |$t_0^+$|⁠, a temperature |$T_H^+$| is applied all along the base of the model: both the viscosity and density evolve and the system progressively destabilizes. The choice of |$T_H^+$| depends on the choice of geodynamic setting that one may reasonably assume: two extreme cases can be considered, either a moderately thickened orogenic crust (H = 45 km) over a thin lithospheric mantle, thus |$T_H^+=1000\,^{\circ}{\rm C}$| , or an orogenic crust that has been overthickened for at least 20 Myr and has thus partially ‘relaxed’ thermally, with H = 60 km and |$T_H^+=900\,^{\circ}{\rm C}$|  (England & Thompson 1986). Since a more detailed parametric study will be carried out in a forthcoming contribution, we only present here a representative case together with a small selection of complementary tests (Table 3). This representative case has H = 60 km and |$T_H^+=1000\,^{\circ}{\rm C}$| (case 1H₅₀), and produces a reasonably good fit with Naxos’s timing and size of convective cells.

Table 3.

Test cases for RB convection in Naxos. Model thickness H, viscosities μ_{bottom; top} and densities |$\rho ^0_{1}$| at 300 °C. The cycles, timing and dimensions are those produced at the end of the model run.

Cases and height	Viscosities	Densities	|$T_H^+$|	Cycles period	Cycles size
(km)	(Pa.s)	(kg m–3)	(°C)	(Ma)	( x × z km)
1H50	|$\mu _{10^{17};10^{21}}$|	2700	1000	1–2	7 × 15 to 15 × 37
2H35	|$\mu _{10^{17};10^{21}}$|	2700	1000	0.7–2	3 × 6 to 12.5 × 20
3H35	|$\mu _{10^{17};10^{21}}$|	2800	1000	0.8–2	10–24
4H50	|$\mu _{10^{18};10^{21}}$|	2700	900		No convection

With |$T_H^+$| applied at the base of the model domain, we determine a Rayleigh number based on the properties of the fluid at average temperature 600 °C. With |$T_H^+=1000\,^{\circ}{\rm C}$|⁠, Ra = 5 × 10⁴. The Prandtl number is then estimated based on the lowest viscosity of the system, 10¹⁷ Pa.s: Pr = 3.7 × 10⁴. This value is lower than realistic geodynamic systems but Krishnamurti (1970) showed that when Pr > 100, the dynamics of the system only depends on the Rayleigh number. Within the Boussinesq approximation, densities span the range 2645−2700 kg m⁻³ after 12.7 Myr.

Results are presented in Fig. 13: the vertical profiles of temperature, density and viscosity (Figs 13a–c), the temporal evolution of some particles (Fig. 13d), and the flow pattern (Fig. 13e). Various sizes of convection cycles appear depending on the initial position of particles. Variations in density and viscosity follow variations in particle temperature. The general flow regime displays 2-D rolls, as predicted by Krishnamurti (1970), with cell sizes of the order of 10 × 20 km and cycle periods ranging from 1 to 2 Myr. These numbers are in the range of those estimated by Vanderhaeghe et al. (2018).

Three other test cases (Table 3) show the following: for a greater viscosity 10¹⁸ Pa.s at H = −60 km depth (case 4H₅₀), the Rayleigh number becomes too low for convection to develop. Decreasing the minimum viscosity below 10¹⁷ Pa.s generates a more vigorous convection, but the convective periods reduce and cell sizes increase, not fitting Naxos’s data anymore. Decreasing the thickness of the model domain decreases the sizes of the convective cells (case 2H₃₅). Finally, increasing the reference density by 100 kg m^–3 only slightly increases cell sizes (case 3H₃₅).

6.2 Naxos Episode 2, diapirism

For the diapirism episode, we consider that the middle and lower crust behave as a three layers system like that of Weinberg (1992). Temperature is not considered here, but from the episode 1 presented above, we choose the top of the model domain to lie 20 km below the top surface, where it becomes significantly less viscous. Underneath, the model domain extends down to –45 km depth, and is made of three layers similar to Weinberg & Schmeling (1992)’s case II: ρ₁ = ρ₂/0.9 = ρ₃/0.89 and μ₂ = μ₃ = μ₁/33, h₂ = h₃ = h₁/8, with an Archimedes number Ar = 0.3 (see setup in Fig. 12b, with parameters given in the caption). The corresponding values stand within the range of variations of Earth’s partially molten middle and lower continental crust (e.g. Rosenberg & Handy 2005; Vanderhaeghe 2009); Hacker et al. 2015).

Figure 12.

Naxos setups: (a) convection and (b) diapirism episodes.

Figure 13.

Naxos RB episode 1, OpenFOAM simulation on a 125 × 125 × 87 grid (see parameters in main text). (a,b,c) vertical profiles of laterally averaged temperature, density and viscosity. (d) Temporal evolution of the depth of 4 passive markers also plotted in e). (e) 3-D view of the isotherms after 12.9 Myr. Convection occurs over a height ∼30 km, producing 4–5 cycles within 7 Myr. The black markers are located along the vertical section z= 16 km, purple and red markers along z= 20 km and blue markers along z= 26 km.

Figs 14(a), (c) and (e) show a vertical section of the system, Figs 14(b), (d) and (f) show a 3-D view of the system, in which ‘domes-in-domes’ structures develop. One of the large grey diapirs which develop within layer 2 is about 10 km wide while white ‘blobs’ inside it are about 2 km wide (see Fig. 14c). The wavelength is about 12.5 km for the first destabilization event (4 diapirs) and 16 km for the second destabilization event (three diapirs). Both little and large diapirs (grey and white layers 2 and 3) rise imbricated, by about 15 km within 0.5 Myr, which corresponds to a velocity of the order of 3 cm yr^–1.

Figure 14.

Naxos RT episode 2, simulation made with OpenFOAM on a 200 × 200 × 140 grid. Total thickness H = 35 km, horizontal widths L = 50 km, reference viscosity μ₁ = 3.3 × 10¹⁸ Pa.s, μ₂ = μ₃ = μ₁/33, densities ρ₁ = ρ₂/0.9 = ρ₃/0.89=2537 kg m^–3 and heights h₂ = h₃ = h₁/8. (a, c, e) display the vertical section along y = 25 km, (b, d, f) display the 3-D interfaces, at three time steps. Diapirs develop within about 1 Myr and reach the upper model boundary within another ≃0.5 Myr.

The modelled widest dome sizes are in agreement with L_Na, the size of the large domes estimated by Vanderhaeghe et al. (2018), and the smallest sizes are similar to d_Na, that of the estimated small domes. The growth rates are also found equivalent to those estimated for Naxos, allowing us to conclude that this model configuration is compatible with the physical state of Naxos’ partially molten crust at some point around 16 Myr.

6.3 Naxos application, discussion for future work

Our preliminary models provide a quantitative description (convection velocities and dome sizes) of the development of gravity instabilities in the orogenic crust of Naxos around 16 Myr, in accordance with the estimates by Vanderhaeghe et al. (2018). These models already yield some insight:

• The relatively low viscosities of the middle-to-lower crust that lead to gravitational instabilities (1 − 5 × 10¹⁷ Pa.s below 20 km depth), given our choices of initial geotherm, densities and basal temperature, imply that this crustal domain behaved according to a dominantly wet quartz composition. Such viscosities match those obtained with rock experiments on hydrated-quartz at mid-crustal conditions (Bürgmann & Dresen 2008). In fact in Naxos Island, the migmatite dome is mantled by metapelites, micaschists and marble layers that are made of quartz, micas and calcite displaying relatively soft rheological behaviour comparable to wet quartz (Kruckenberg et al. 2011).

• The modelled dynamics occurs at a deformation rate of the order of 10⁻¹³−10⁻¹⁴s⁻¹ corresponding to several cm per year of vertical mass transfer. This sets a perspective on the relative influence of far field tectonic boundary conditions, progressively switching from compressional to extensional over a time period of about 10 Myr throughout the Cyclades orogeny. Such a transition in tectonic kinematics is also invoked in other orogenic contexts, and raises the question of how heat from the mantle below can be maintained over equivalent time periods (England & Thompson 1986; Ueda et al. 2012; Gerbault et al. 2018). This subject needs further work beyond the scope of the present contribution.

• The convection model (episode 1) shows that several cycles of convection of distinct timing occur, depending on whether the edges or the core of the rising crustal plumes are being sampled. These sampled locations cover distinct ranges of temperature: the large cycles cover a range of 670−920 °C (purple marker in Figs 13d–g), whereas the small cycles remain at 810 ± 10 °C (orange marker in Figs 13d–g). Since zircons typically crystallize and dissolve as they progressively pass through the ∼[700−900] °C range (Guergouz et al. (2018) and references therein), both small and large cycles may then explain the pattern that the zircons of Naxos appear to have gone through. Complementary data on sampled zircons across the domes of Naxos will help to better constrain their dynamics.

Further work is required to understand better how Naxos’s crustal system would have progressively switched from convection to polydiapirism, or vice versa. Although this switch may have been linked to a change in kinematic or thermal boundary conditions, more complex rheologies accounting for melting or internal heating were also likely at play. Some of these processes will be adressed in a forthcoming paper. One may also add that sub-scale two-phase flow of low viscosity and low density melt fluids should also be accounted for. Nonetheless, the overall mass balance associated with the evolution of this crustal domain would also require to make hypotheses on (i) the amount of eroded material from the top surface and (ii) the amount of material that has sunk down below Moho levels. These issues are clearly out of the scope of this study, and we refer to other studies that discuss them (Weinberg 1997; Jull & Kelemen 2001; Burg & Vigneresse 2002; Gerya et al. 2008; Gerbault et al. 2018; Riel et al. 2016; Cao et al. 2016; Piccolo et al. 2019; Schmeling et al. 2019).

7 CONCLUSIONS

We have tested two VOF methods without interface reconstruction, for four well documented Rayleigh–Taylor and Rayleigh–Bénard instabilities systems, and the results are shown to be in good agreement with analytical solutions and reference numerical and experimental results. Both JADIM and OpenFOAM codes were shown to be able to reproduce the various convection regimes that are predicted in a temperature-dependent two-layer system heated from below, for a reasonable range of parameters, both in two and three dimensions. To our knowledge such detailed comparisons of VOF implementations have not been produced in the literature before. Using either a harmonic or an arithmetic viscosity averaging at interfaces leads to different structures in the deformed layers, even though the models match theoretical predictions for upwelling velocities and growth rates. The compositional fields tracking method used here via the geosciences code ASPECT appears more sensitive to initial perturbations of interfaces than the two VOF codes. The VOF codes exhibit less diffusion at interfaces and mass is better conserved in comparison to ASPECT (by several orders of magnitude).

In terms of codes performances, OpenFOAM is about two times faster than JADIM. In 2-D, the compositional field tracking method available in ASPECT is about two times faster than OpenFOAM, mainly because the viscous term is solved with an implicit method while it is solved semi-implicitly in both VOF methods. Despite a better parallel efficiency of ASPECT, OpenFOAM appears faster than ASPECT in the specific 3-D Rayleigh–Bénard case that was tested (given the default computational options that we used).

Here, the VOF method has been applied for the first time to 3-D three-layer Rayleigh–Taylor systems of crustal geodynamics. While some properties like the initial growth rates are similar between 3-D and 2-D models, differences such as the timing of destabilization reach about 20 per cent. Therefore, we conclude that 3-D simulations are necessary in order to study precisely the structural evolution of gravitational instabilities of stratified crustal layers, in comparison to field data.

A preliminary application to the gneiss dome and subdomes exposed in Naxos is addressed by testing independently a convection setup and a diapirism setup. Both mechanisms or episodes, modelled with OpenFOAM, show good agreement with previous analytical estimates based on zircon geochronology and structural field geology: both mechanisms are likely to have been at play in the formation of the polydiapiric structures observed on Naxos Island. This provides insight into the geodynamic setting and rheological properties of this hot crustal domain over a period of about 15 Myr. However, several questions remain about the actual transition or combination of both these episodes. A follow-up study aims at investigating how both episodes could operate from a common starting configuration, as well as testing other first order parameters that would have controlled the evolution of this hot crustal domain.

To summarize, the present work shows that the VOF method is a promising tool for studying crustal diapirism and convection, aiming at a better understanding of the thermomechanical processes responsible for the exhumation of partially molten orogenic crust.

SUPPORTING INFORMATION

ASPECT_SETUP.zip

Fig6a_3D_CaseII_JADIM.mp4

Fig7a_3D_Case1_Weinberg_OpenFOAM.mp4

Fig10b_ASPECT.mp4

Fig10b_JADIM.mp4

Fig10b_OpenFOAM.mp4

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

This work benefited from fruitful discussions with Anne Davaille, Yves D’Angelo and Annaïg Pedrono. An anonymous reviewer, René Gassmöller and editor Gael Choblet are warmly thanked for constructive reviews. Some of the computational time was provided by the super computing centre CALMIP (Toulouse, project P152). We are grateful to the Institut National Polytechnique de Toulouse (INPT) which supported this work. This paper is part of the PhD thesis of the first author. The authors would like to thank Rached Abdesslem and Liwei Deng who contributed to the initiation of this collaborative project through their master’s degree internships. We thank the Computational Infrastructure for Geodynamics (https://geodynamics.org/) which is funded by the National Science Foundation under award EAR-0949446 and EAR-1550901 for supporting the development of ASPECT. Funding from INSU for a 3 days field mission to Naxos is acknowledged.

APPENDIX A: WAVELENGTH AND GROWTH RATE PREDICTIONS FOR TWO- AND THREE-LAYER SYSTEMS

A1 Ramberg (1981a) analytical predictions

The vertical position y_i of a layer interface i is assumed to evolve as |$y_i(t) = A_i e^{\kappa _i q_i t}$|⁠. Ramberg (1981a) derived the following linear system for a two or three layers system to solve for vertical displacements y₁, y₂:

$$\begin{eqnarray*} \quad {\begin{bmatrix}c_{11} & \quad d_{12} & \quad e_{13} & \quad f_{14} \\ i_{21} & \quad j_{22} & \quad k_{23} & \quad l_{24} \\ a_{31} & \quad b_{32} & \quad c_{33} & \quad d_{34} \\ g_{41} & \quad h_{42}& \quad i_{43} & \quad j_{44} \\ \end{bmatrix}} {\begin{bmatrix}v_1 \\ u_1 \\ v_2 \\ u_2 \\ \end{bmatrix}} = {\begin{bmatrix}0 \\ q_1y_1 \\ 0 \\ q_2y_2 \\ \end{bmatrix}} . \end{eqnarray*}$$

(A1)

All coefficients c₁₁ to j₄₄ depend on ϕ_i, μ_i and boundary conditions. v_i and u_i are the velocities of layer i in horizontal and vertical directions, respectively. q_i is defined as |$\displaystyle {q_i = \frac{(\rho _i - \rho _{i+1})g h_{i+1}}{2 \mu _{i+1}}}.$|

For example, for a two-layer system, the dimensionless growth rate κ is:

$$\begin{eqnarray*} \displaystyle { \kappa = \frac{-d_{12}}{c_{11}j_{22}-d_{12}i_{21}}} \end{eqnarray*}$$

(A2)

with:

$$\begin{eqnarray*} \displaystyle { c_{11} = \frac{\mu _12\phi _1^2}{\mu _2(\cosh (2\phi _1-1-2\phi _1^2)}-\frac{2\phi _2^2}{\cosh (2\phi _2)-1-2\phi _2^2} } \end{eqnarray*}$$

(A3)

$$\begin{eqnarray*} \displaystyle { d_{12} = \frac{\mu _1 (\sinh (2\phi _1)-2\phi _1}{\mu _2(\cosh (2\phi _1)-1-2\phi _1^2)}+ \frac{\sinh (2\phi _2)-2\phi _2}{\cosh (2\phi _2)-1-2\phi _2^2}} \end{eqnarray*}$$

(A4)

$$\begin{eqnarray*} \displaystyle { i_{21}=\frac{\mu _1 \phi _2(\sinh (2\phi _1)+2\phi _2)}{\mu _2(\cosh (2\phi _1)-1-2\phi _1^2)} + \frac{\phi _2(\sinh (2\phi _2)+2\phi _2)}{\cosh (2\phi _2)-1-2\phi _2^2} } \end{eqnarray*}$$

(A5)

$$\begin{eqnarray*} \displaystyle { j_{22} = \frac{\mu _1 2 \phi _1^2 \phi _2}{\mu _2(\cosh (2\phi _1)-1-2\phi _1^2)} - \frac{2\phi _2^3}{\cosh (2\phi _2)-1-2\phi _2^2}} \end{eqnarray*}$$

(A6)

and

$$\begin{eqnarray*} \phi _i = \frac{2 \pi h_i}{\lambda }. \end{eqnarray*}$$

(A7)

A2 Burg et al. (2004) analytical predictions

Burg et al. (2004) provided relations to find the dominant wavelength λ_dom and growth rate κ for a semi-infinite two-layer system:

$$\begin{eqnarray*} \displaystyle {\frac{{\lambda _{dom}}}{h_2}=0.79 \left(\frac{\mu _t}{\mu _b} \right)^{1/3}}, \end{eqnarray*}$$

(A8)

$$\begin{eqnarray*} \displaystyle {\kappa = \frac{\sigma }{q}=0.53 \times 2 \times \left(\frac{\mu _t}{\mu _b} \right)^{-2/3}}, \end{eqnarray*}$$

(A9)

where h₂ is the depth of the bottom layer, and μ_b and μ_t are the viscosity of the bottom and top layers, respectively.

APPENDIX B: DIMENSIONLESS SYSTEM OF EQUATIONS

All equations in the codes are solved in dimensional form, but they are formulated here after in their dimensionless form and expressed in terms of the Archimedes number (Appendix  C):

$$\begin{eqnarray*} \frac{\partial C_i}{\partial t^* } + {\bf U}^* \cdot \nabla C_i = - \nabla \cdot ( {\bf U}_r^* C_r ), \end{eqnarray*}$$

(B1)

$$\begin{eqnarray*} \nabla \cdot {\bf U}^* = 0, \end{eqnarray*}$$

(B2)

$$\begin{eqnarray*} \frac{\partial {\bf U}^* }{\partial t^* } + {\bf U}^* \cdot \nabla {\bf U}^* &=& - \nabla P^* + {\bf g}^* \nonumber \\ &&+ \nabla \cdot \left[ \frac{2}{Ar} (\nabla {\bf U}^* + (\nabla {\bf U}^* )^T)\right], \end{eqnarray*}$$

(B3)

$$\begin{eqnarray*} \frac{\partial T^* }{\partial t^* } + {\bf U}^* \cdot \nabla T^* = \nabla \cdot \left(\frac{2}{Ra} \nabla T^* \right), \end{eqnarray*}$$

(B4)

with: U* = U/qh; P* = P/ρ(qh)²; g* = g/(q²h).

APPENDIX C: INFLUENCE OF ARCHIMEDE’S NUMBER

Fig. C1 shows the dimensionless mean velocity of the two-layer Rayleigh–Taylor system depicted in Fig. 1a) (with γ = 100 and h₂/H = 0.2) for Ar ≥ 1 (a) and Ar ≤ 1 (b) as a function of dimensionless time τ. For Ar ≥ 1, the dimensionless mean velocity varies according to Ar. For Ar ≤ 1, dimensionless mean velocities are equal and follow the linear stability theory.

Figure C1.

Time evolution of the mean velocity in the two layer system for various Archimedes numbers. Velocity is scaled by |$\mathcal {U}$| and time is scaled by |$\mathcal {T}$| (eq. 11). For comparison, the result of the linear stability theory of Ramberg (1981a) is also plotted (denoted by LST).

However, at small τ (<10), the velocity differs due to the initiation of destabilization. Thanks to these results, we expect that the dynamics of a geological system (Ar typically around 10⁻²⁰) are equivalent to those of a system with Ar = 1. We can then use a viscosity around 10³ Pa.s instead of 10²⁰ Pa.s, which accelerates our calculations. This result was shown in a Rayleigh–Taylor case but it is likely to be valid for any configuration as long as the local instantaneous Ar ≤ 1. If Ar ≥ 1, inertial effects may become significant and change the flow dynamics. Note that in the present context, the Archimedes number is nothing but a Reynolds number using a velocity scale which may be viewed as a Stokes velocity. Therefore, it is reasonable to say that having Ar ≤ 1 locally leads to negligible inertia effects.

APPENDIX D: VAN KEKEN’S TEST: VISCOSITY AVERAGING

The main text displays van Keken’s bench for viscosity contrasts γ = 100 (Fig. 1) and an arithmetic viscosity averaging. The secondary instability on the top right is the most delicate to reproduce. Here we present the results of JADIM, OpenFOAM at γ = 1 and 10 with arithmetic viscosity averaging (Fig. D1). In addition, γ = 10, 100 with an harmonic law are displayed Fig. D2, with the addition of ASPECT results.

Figure D1.

Evolution of van Keken’s system for viscosity ratios (a) γ = 1 thus no averaging, (b) γ = 10 and arithmetic averaging, at |$t/\mathcal {T}=50$|⁠, 100, 150. Van Keken's snapshots are modified from Figs. 2 and 4 of van Keken et al. (1997). Resolutions 91×100 for JADIM and OpenFOAM.

Figure D2.

Evolution of the two-phase RT system for a viscosity ratio γ = 100 at times |$t/\mathcal {T}=50$|⁠, 100, 150. From top to bottom, results from van Keken et al. (1997) (Figs. 4 and 6), JADIM and OpenFOAM with harmonic viscosity averaging at resolution 91 × 100. Bottom-right displays ASPECT results, with harmonic viscosity and resolution 91 × 100.

Figure D3.

Comparison of the time evolution of the mean velocity in the Rayleigh–Taylor two-layer system by van Keken et al. (1997), with JADIM and OpenFOAM, with harmonic viscosity averaging and a resolution 100 × 100. Velocity is scaled by |$\mathcal {U}$| and time is scaled by |$\mathcal {T}$|⁠, cf. eq. (11). |$\gamma = \frac{\mu _1}{\mu _2}$|⁠.

APPENDIX E: INFLUENCE OF MESH RESOLUTION, VISCOSITY INTERPOLATION AND INTERFACE PERTURBATION ON RAYLEIGH–TAYLOR INSTABILITIES

We plot in Fig. E1(b) the time evolution of the velocity of the maximum vertical location of interface 1–2 with or without perturbing the two interfaces with a random perturbation in JADIM. For t/τ ≤ 180, the growth rate of the unperturbed interface is higher than that of the perturbed interface but the velocity value is lower. For t/τ ≥ 300, the growth rates become similar. The shape of interface 2–3 however, is significantly different depending on the type of initial perturbation (see snapshots in Fig. E1). Without initial perturbation both instabilities have a roughly similar wavelength and hence interface 2–3 deforms as an unique diapir. This is different from the prediction of the linear stability theory (Fig. 5d). Alternatively, the results for the simulation with perturbation are in agreement with the theory (Figs 5c and d). It turns out that in this case interface 2–3 has a wavelength very sensitive to the initial perturbation of interfaces. This may be explained by how close the growth rates of the two interfaces are. This is not the case for cases I (Fig. E1a) and III (Fig. E1c) where the shape and dynamics of the interfaces is roughly independent of the initial perturbation of the interfaces.

Figure E1.

Case I (a), II (b) and III (c) of Weinberg & Schmeling (1992): time evolution of the velocity of the maximum height of interface 1–2: with or without an initial random perturbation at interfaces (JADIM code, harmonic viscosity averaging), compared to linear stability theory (Ramberg 1981a). Time is scaled by |$\mathcal {T}$| and |$t/\mathcal {T}=0$| is set as the time when one of the interface has raised by a distance of 3 × 10⁻⁵H.

In the main text we display the results of Rayleigh–Taylor three-phase systems with harmonic viscosity interpolation. Below, in Figs E2 and E3 we show complementary tests comparing velocities between JADIM, OpenFOAM and ASPECT for various mesh resolutions and viscosity averaging.

Figure E2.

Influence of mesh resolution on Weinberg & Schmeling (1992) case II: time evolution of the mean velocity for (a) JADIM, (b) OpenFOAM and (c) ASPECT, with harmonic viscosity averaging, γ = 100. Note that an initial perturbation at the interfaces was imposed for JADIM and ASPECT, not for OpenFOAM.

Figure E3.

Three-phase system by Weinberg & Schmeling (1992), with an arithmetic or a harmonic interpolation of viscosity at interfaces, with JADIM: Case II, the top layer is the densest, γ = 100. Time is scaled by |$\mathcal {T}$| and |$t/\mathcal {T}=0$| is set as the time when one of the interface has raised by a distance of 3 × 10⁻⁵H.

The time evolution of the maximum vertical velocity at interface 1–2 using an arithmetic or a harmonic viscosity averaging are compared in Fig. E3. For 0 < t/τ < 200, the arithmetic averaging exhibits a higher growth rate than the harmonic averaging, but as t/τ > 250, the growth rates become similar.

The evolution of layers 1 and 3 is roughly similar independently of the viscosity averaging. Note however, that at time t/τ = 497 with arithmetic averaging, the diapir of layer 2 is slightly smaller than that with harmonic averaging. Indeed, the harmonic averaging favours a behaviour controlled by the lowest viscosity at an interface, thus it is easier for layer 2 to invade layer 1.

Table E1.

Relative error of the dimensionless growth rate κ as a function of mesh resolution. The case of resolution 240 × 100 is arbitrarily taken as reference for JADIM and OpenFOAM.

	Viscosity averaging	harmonic	harmonic	harmonic	harmonic	arithmetic
	Grid size	60 × 25	90 × 37	120 × 50	240 × 100	240 × 100
JADIM	1st instability	0.1  %	5  %	3  %	REF	1  %
	2nd instability	25  %	2  %	1.2  %	REF	2  %
	Grid size	60 × 25	90 × 37	120 × 50	240 × 100	240 × 100
OpenFOAM	1st instability	5  %	5  %	5  %	REF	0.8  %
	2nd instability	22 %	4.5  %	4.5  %	REF	0  %

APPENDIX F: MESH CONVERGENCE TESTS WITH JADIM FOR Vatteville et al. (2009)’S CASE

Figure F1.

Rayleigh–Bénard problem (Vatteville et al. 2009): comparison of velocity between laboratory (Vatt-Exp) and numerical experiments with JADIM, testing the influence of grid size. Velocity is scaled with V_Stokes = αgΔTH²/ν_max, time is scaled with τ = ν_max/(αgΔTH), and distance with the domain height H.

APPENDIX G: MASS ERROR FOR DIFFERENT SETUP CASES

Figure G1.

Relative mass error ΔM as a function of time for JADIM, OpenFOAM and ASPECT: (a) Two-layer Rayleigh–Taylor system (γ = 100, harmonic viscosity averaging, resolution 91 × 100); (b) Three-layer 2-D Rayleigh–Taylor system (solid line: layer 3; dotted line: layer 1 (2 for OpenFOAM); case II, harmonic viscosity averaging, resolution 240 × 100); (c) Two-layer Rayleigh–Bénard system (harmonic viscosity averaging, resolution 90 × 90 × 44); (d) Same as (b) in 3-D (resolution 120 × 100 × 120).

APPENDIX H: CODES PERFORMANCES

Figure H1.

Strong (a) and weak (b) scaling for the two-layer Rayleigh–Bénard problem of Le Bars & Davaille (2004), case II, for JADIM, OpenFOAM and ASPECT with 10 time-steps. Dotted lines indicate optimal parallelization. For the strong scaling (a) a same 120 x 120 x 60 grid model is taken for all the simulations. For the weak scaling (b), the CPU time is scaled over the CPU time spent by a same basic 30x30x30 model run on 1 core.