Using the level set method in geodynamical modeling of multi-material flows and Earth ’ s free surface

The level set method allows for tracking material surfaces in 2-D and 3-D flow modeling and is well suited for applications of multi-material flow modeling. The level set method utilizes smooth level set functions to define material interfaces, which makes the method stable and free of oscillations that are typically observed in case step-like functions parameterize interfaces. By design the level set function is a signed distance function and gives for each point in the domain the exact distance to the interface as well as on which side it is located. In this paper we present four benchmarks which show the validity, accuracy and simplicity of using the level set method for multi-material flow modeling. The benchmarks are simplified setups of dynamical geophysical processes such as the Rayleigh–Taylor instability, post-glacial rebound, subduction and slab detachment. We also demonstrate the benefit of using the level set method for modeling a free surface with the sticky air approach. Our results show that the level set method allows for accurate material flow modeling and that the combination with the sticky air approach works well in mimicking Earth’s free surface. Since the level set method tracks material interfaces instead of materials themselves, it has the advantage that the location of these interfaces is accurately known and that it represents a viable alternative to the more commonly used tracer method.


Introduction
Accurate modeling of geodynamical processes involving large deformation, e.g., mantle flow, subduction evolution or slab tearing, is a key research goal in computational geody-namics.Since the early simplified two-dimensional isothermal model configurations (Gurnis and Hager, 1988;Christensen, 1996), model complexity and especially the number of materials present in numerical models has dramatically increased as seen in recent three-dimensional thermomechanically coupled models that include multiple materials and phase changes as well as surface deformation and complex rheologies (e.g., van Hunen and Allen, 2011;Duretz et al., 2014).For instance multiple materials are important to investigate the influence of an oceanic crust on the decoupling of subducting and overriding plates as well as on the buoyancy of the subducting slab (Běhounková and Čižková, 2008;van Hunen and van den Berg, 2008;Androvičová et al., 2013); they are also important for research involving subduction termination by continental collision (e.g., Baumann et al., 2010;Magni et al., 2012).Other studies focus on the influence of complex rheologies on slab dynamics (Billen and Hirth, 2007;Andrews and Billen, 2009); finally, the aim of many recent studies has been to investigate the influence and response of a free surface in subduction modeling (e.g., Schmeling et al., 2008;Gerya et al., 2009;Quinquis et al., 2011;Duretz et al., 2011).
The models described above invariably require the ability to track different materials and material interfaces throughout the model domain.Within the community several different methods, based on either Lagrangian or Eulerian modeling frameworks, are used.In Lagrangian finite element codes, the mesh is deformable and element boundaries are often aligned with the material interfaces.The elements thus track the materials through the model.However, when large deformation is modeled, e.g., when following subduction evolution, the B. Hillebrand et al.: Geodynamical modeling using the level set method mesh may become too distorted and remeshing is required.This is computationally expensive, results in unwanted numerical diffusion and constitutes an important drawback for models with large deformation.
In Eulerian codes the mesh is fixed.Because of this elements do not track materials and a method for material tracking is needed.The two more commonly used methods in computational geodynamics are the marker-in-cell technique and the phase field method.
Tracers (particles, markers) are widely used in the geodynamical community (e.g., Tackley and King, 2003;van Hunen and Allen, 2011;Duretz et al., 2012).These Lagrangian particles are advected with the flow and carry material properties such as density and viscosity.Velocity equations are solved on the finite element mesh and the velocities used to advect the particles are obtained from interpolation.One element generally contains several up to 100s of particles.This method easily allows for advection of multimaterial fields using an Eulerian mesh and is potentially nondiffusive (Tackley and King, 2003).While the tracer method tracks materials and has as an advantage that it can record its history, it does however not track the interface between the materials.The interface position remains approximate and is known with an uncertainty on the order of local tracer distance.However, the tracer method tracks materials, it does not track the interface between the materials.The interface position remains approximate and is known with an uncertainty on the order of local tracer distance.Furthermore, the tracer method becomes increasingly expensive in 3-D.For instance, in the 2-D models of Crameri et al. (2012) the different codes used between ten and hundreds of particles per element/cell.In 3-D this translates to several dozens up to thousands of particles per element/cell, i.e., possibly billions in total in the case of large 3-D simulations.The computational as well as the memory costs would then become huge, requiring that the code is highly parallel and scales up to hundreds of cores or more.
In the phase field method (Lenardic and Kaula, 1993;Van Keken et al., 1997;Kronbichler et al., 2012), materials are assigned a number, and the composition of the fluid at a given node of the grid is given by a field containing the various fractions of the different material components.This field is then advected using a stabilized advection equation.The phase field vector is only defined on the nodal points of the mesh; thus, the computational costs increase proportionally to the increase in nodal points.However, such a phase field will contain sharp contrasts between the different phases within elements and the advection of the phase field requires complex stabilization schemes (Lenardic and Kaula, 1993).
In this paper we explore a third method, the level set method that, instead of tracking materials, is geared to track the material interfaces.The method is based on defining a level set function (generally signed distance) which is zero at the target interface and positive on one side and negative on the other side.This signed property is used to identify the dif-ferent materials.Similar to the phase field method the level set function is defined on the nodal points of the elements and the computational costs increase proportional to the increase in nodal points.In contrast with the phase field method, the level set method does not involve step-like discontinuities but instead represents fields with a smooth (signed distance) function (the level set function).The level set method provides a sharp location for the interface.
The level set method has not often been used in the geodynamical community.Notable exceptions are Hale et al. (2007), Gross et al. (2007) and Bourgouin et al. (2007) which are focussed on lava dome growth and/or mantle plumes, Suckale et al. (2010) on bubbles dynamics in volcanic conduits, Zlotnik et al. (2008) on gravitational instabilities and Hale et al. (2010) on slab tear faults.In Braun et al. (2008) a level set method is presented which is based on a 3-D set of triangulated points, which makes it a hybrid between tracers and level set functions.The level set method is primarily used in other fields of computational science such as twophase flows (Oka and Ishii, 1999) and fluid dynamics (Rao et al., 2011).An overview of the method and applications can be found in Osher and Fedkiw (2001).
In this paper we present four benchmarks of increasing complexity and end with applications to modeling of subduction and slab detachment.Two of the presented benchmarks include deformation of Earth's free surface.In Eulerianbased codes this is generally modeled either by ALE (arbitrary Lagrangian-Eulerian; Fullsack, 1995;Thieulot, 2011) methods or the so-called sticky air approach (Schmeling et al., 2008;Crameri et al., 2012).In ALE the top layers of elements can deform vertically.The sticky air approach entails that an air layer of low viscosity and zero density is put atop the surface.This causes the Earth's surface to become a boundary between two materials inside the model domain which we track using the level set method.
The purpose of our paper is to demonstrate the use of the level set method in various geodynamic applications and to demonstrate the applicability of the presented approach to more complex geodynamical processes.

Methods
All experiments in this paper are mechanical models internally driven by density perturbations.They do not include any temperature effects.We use the finite element modeling package SEPRAN (Segal and Praagman, 2005) and solve for mass conservation of an incompressible fluid: and the Stokes equation describing a force balance: Here v is velocity, P dynamic pressure, ρ density and σ d the deviatoric stress tensor.In all benchmarks the density is a function of the material that is advected through the domain which is tracked by means of the level set method.The equations are solved on linear triangular elements.The modeling package SEPRAN has been used in geodynamical modeling for many years(e.g., Čížková et al., 2007;van Hunen and van den Berg, 2008;Chertova et al., 2012;Androvičová et al., 2013).

Introduction
The level set method is a well-researched interface tracking technique which was originally devised by Osher and Sethian (1988).It tracks an interface by defining it as the zero valued isocontour of a smooth function.In the last two decades several improvements and variations have been presented by several authors such as reinitialization (see below) extension velocities (e.g., Adalsteinsson and Sethian, 1999;Chopp, 2009), local level set methods (e.g., Sethian, 2001), hybrid particle level set methods (e.g., Enright et al., 2002;Samuel and Evonuk, 2010), variational level set method (e.g., Duan et al., 2008) and the level set method combined with volume of fluid method (e.g., Fedkiw et al., 1999;Pijl et al., 2008).
The level set method implementation presented here is global and uses reinitialization to keep the level set function smooth.If denotes the interface that is to be associated and tracked with a level set function φ, and is a bounded region, bounded by just the interface or the interface and the boundaries of the model domain, φ will be defined as (Osher and Fedkiw, 2001): The level set function is advected by means of the advection equation: This equation is solved using a Crank-Nicolson integration scheme in combination with the streamline upwind Petrov-Galerkin (SUPG) upwind scheme (Brooks and Hughes, 1982).The level set function is generally chosen to be a signed distance function which means that |∇φ| = 1 everywhere.The function value indicates on which side of the interface a point is located (negative or positive) and this is used to identify materials.Because the level set function is a signed distance function, its value is also the distance to the interface.

Reinitialization
The level set function is advected with a velocity field resulting from the buoyancy forces.This velocity field does not necessarily preserve the signed distance quality of the level set function.However, it has been shown by several authors (Sussman et al., 1995;Min, 2010) that it is important for the level set function to stay smooth in the vicinity of the zero level set (at least Lipschitz continuous; Osher and Fedkiw, 2001).Therefore, the level set function is corrected so that it remains a smooth function without moving the zero isocontour and thus the interface itself.Sussman et al. (1995) introduced a method called reinitialization which exploits the signed distance quality of the level set function.The reinitialization process involves solving the following equation: This equation specifies a correction for the value of φ if |∇φ| = 1.∂τ is a pseudo time step and sign(φ) is the one dimensional signum (or sign) function, and 0 at φ = 0, for which generally a smooth approximation is used.Equation (5) does not need to be solved every time step.We use an error criterion to determine whether reinitialization is needed and then solve several reinitialization iterations until our criterion is satisfied.The error is calculated by taking the average of the deviation of the absolute gradient of φ from one of all the nodal points.The number of iterations needed depends on the choice of ∂τ , sign(φ) and the choice of the error criterion.For the smoothened sign function we use The value of constant C is arbitrarily chosen.A high value results in a slower (more iterations) but more stable reinitialization process while a low value has the opposite effect.In our models C = 15 proved to be a practical compromise between speed and stability.The determination of ∇φ is important for the reinitialization procedure and it needs to be robust for both small and large-scale variations.We therefore use a second order ENO (essentially non-oscillatory) scheme for the space derivative (Osher and Shu, 1991;Jiang and Peng, 2000).We also use a third order TVD (total variation diminishing) Runge-Kutta scheme for the pseudo time integration of Eq. ( 5) (Gottlieb and Shu, 1998).

Usage of the level set method
We use the level set method to track the interface between two different (geological) materials, but we note that the application can involve any chosen surface.Every interface is described by its own level set function.Since the level set function is defined such that its zero isocontour coincides with the interface between two materials, one material will be where the function is negative and the other where it is positive.For an arbitrary material parameter C this can be written as We coin this the sharp boundary method which results in sharp contrasts of material parameters (density, viscosity, etc.) within an element.One can also introduce a small diffusion zone around the interface (Bourgouin et al., 2006), hereafter called the diffuse boundary method.It is important to note that this does not mean that the location of the interface is no longer known.The location is still exactly known.
Here α = 1 and h represents one element size.When Eq. ( 8) is used to smooth density across the interface, C is simply the density value; however, when viscosity is smoothened across the interface, C is the exponent of viscosity (i.e., smoothing the logarithm of viscosity).The level set function is evaluated at the nodal points of the elements.Density and viscosity are thus assigned to the nodal points and then interpolated to the Gaussian integration points of the finite elements.
A comparison between the two methods is performed with the Rayleigh-Taylor instability benchmark found in Sect. 4.
In the other benchmarks the diffuse boundary method is used.Because of the signed distance quality of the level set function, the zone of the diffuse boundary (2h) follows directly from the level set function values and no additional steps to identify this zone are required.The width (h) of the diffuse boundary is easily changed in case more smoothness is required.

Sticky air approach
Several of the benchmark experiments we conduct will include an approximation of the Earth's free surface using the so-called sticky-air approach (Schmeling et al., 2008;Crameri et al., 2012).This allows modeling of topography changes while using a purely Eulerian code by augmenting the model with a top layer with so-called sticky-air.Since Earth's surface is effectively a stress-free surface this layer of air should exert as little stress on the underlying lithosphere material as possible.Crameri et al. (2012) investigated the viscosity contrast and thickness of the sticky air layer and concluded that for a 100 km thick layer the viscosity of the air layer should be 5 orders of magnitude less than the underlying material.The density of the sticky air layer is set to zero so that it has no pressure effect on the real free surface (the sticky air-lithosphere interface).

Benchmarks
Here we present the model setups of the various benchmarks presented in this paper.All four benchmarks describe multimaterial flow models and some include the modeling of Earth's free surface by means of a sticky air layer.The benchmarks are the Rayleigh-Taylor instability from Van Keken et al. (1997), the post-glacial rebound setup from Crameri et al. (2012), the free subduction benchmark from Schmeling et al. (2008) and the simplified slab detachment setup from Schmalholz (2011).The first benchmark models the overturn of a gravitationally unstable compositional layering and is often used in the geodynamical modeling community.The second and third benchmark focus on the sticky air approach.
The last benchmark setup demonstrates the splitting of one material domain into two.The setups of the benchmarks are illustrated in Fig. 1 and a small description of each follows below.

Rayleigh-Taylor instability
This benchmark represents a buoyancy driven flow and (Fig. 1a) has been performed by several authors with various techniques including tracers (Van Keken et al., 1997;Tackley and King, 2003), level set method (Bourgouin et al., 2006), particle level set method (Samuel and Evonuk, 2010), phase field method (Bangerth and Heister, 2013) and a marker chain method (Van Keken et al., 1997).The benchmark describes an almost square domain of unit height and a width of 0.9142, in which a dense layer overlies a lighter layer.
The interface geometry between the two layers is given by a sine function defined as w(x) = 0.02 cos( πx λ ) + 0.2 with λ = 0.9142.We will compare snapshots at regular intervals with the snapshots of the original article (Van Keken et al., 1997).We will also compare the evolution of the root mean square velocity (v rms ) of the entire domain over time, specifically concentrating on the timing and height of the first peak, which coincides with the rise of the first diapir.The v rms is given by Here V is the volume of the domain.

Post-glacial rebound
This benchmark (Fig. 1c) is used specifically to validate the sticky air approach.It describes a three layer model (air, lithosphere and mantle) with three different viscosities and two different densities (mantle and lithosphere have the same density) and therefore requires two different level set functions at the two-material interfaces.The surface (top of lithosphere) has a prescribed cosine topography with a 7 km amplitude.The system relaxes over time until the topography has reduced to zero.The height of the topography at the left side of the model is measured over time and will be compared to the semi-analytical solution from the original article.

Subduction benchmark
This subduction setup (Fig. 1b) was presented as a benchmark in Schmeling et al. (2008) and this particular setup had been performed by five different codes therein.It involves three different materials: a sticky air layer, an idealized slab which subducted for a 100 km and a mantle.It contains two level set functions which partly overlap: one tracking the interface between the air and the mantle/lithosphere and one tracking the interface between the slab and the air/mantle.The zero level sets of the two level set functions can be seen in Figs. 6 and 7. Due to its negative buoyancy the slab starts to develop rollback and sinks into the mantle.The original paper clearly highlighted difficulties with different choices in tracer-based viscosity averaging schemes due to the entrainment of tracers, a problem we aim to avoid by using the level set method.For comparison we will focus on the depth of the slab tip with time and snapshots through time.

Slab detachment
This setup (Fig. 1d) is from Schmalholz (2011) and is being developed into a community benchmark (Thieulot et al., 2014b).It concerns a two-material model setup comprising a lithosphere with a vertically hanging slab and a mantle.The two materials have different densities and different rheological parameters.Mantle material has a linearly viscous (n = 1, η 0 = 10 21 Pa s) rheology while the slab follows a power-law rheology described by where ˙ is the second invariant of the strain-rate tensor.The following values are adopted: n = 4 and η 0 = 4.75 × 10 11 Pa s.We measure the thickness D of the thinning slab over time.We present our results in the same non-dimensional form as Schmalholz (2011), i.e., nondimensional thickness D d = D D 0 vs. non-dimensional time t d = t t c .D 0 is the initial thickness of the slab (80 km) and t c is the characteristic time.This time is calculated in the following manner: where g is the gravitational acceleration, ρ the density contrast between slab and mantle, H the length of the hanging slab and B = (2η 0 ) −n .The benchmark illustrates the separation of the level set field into two domains.

Rayleigh-Taylor instability
In Fig. 2 the time evolution of both the density and level set field of a 160 × 160 elements model run are shown.This can  be compared with Fig. 2 from Van Keken et al. (1997).All the large-scale features (the two upwellings, the major downwellings) are captured as well as the smaller-scale features such as the small wavelet just behind the front of the first upwelling (Fig. 2c and d).The evolution of the level set field illustrates the signed distance quality of the function.In Fig. 3a the root mean square velocity (v rms ) of the entire domain is plotted vs. time.The first peak corresponds to the first upwelling and the second smaller peak to the second upwelling.
The figure shows the results of four of our models as well as results from the marker chain method from Van Keken et al. (1997).Figure 3c shows a close-up of the first peak.This close-up shows that our results are in good agreement with Van Keken et al. (1997).The first peak is a strong feature across different codes and is therefore examined more precisely in Table 1.It is compared to the marker chain method of Van Keken et al. (1997), the level set method and the particle level set method of Samuel and Evonuk (2010) and the tracer method of Thieulot et al. (2014a).Table 1 illustrates that the first peak occurs earlier with increasing resolution, and that our results resemble the marker chain results of the original paper more than the other published level set method results.However, solely looking at the highestresolution models all presented studies agree on timing and height within 1 %.12).(c) A zoom-in of the first peak in the rms velocity plot of (a).
Figure 3a and c shows the results of a 160 × 160 sharp boundary method run and of a diffuse boundary method run.Although the overall difference is small the diffuse boundary method has a beneficial smoothening effect (Fig. 3c).As previously stated the signed distance quality of the level set function makes such a diffuse boundary method cheap and simple.
This benchmark involves flow of an incompressible fluid.Given the large deformation of the layers, it is appropriate to test the mass conservation of the system.To this end we calculated the relative mass difference defined in Eq. ( 12) and plotted the results in Fig. 3b.
where M 0 is the mass of the entire model domain at t = 0.Over time the models exhibit a variation in total mass.At t = 2000 the relative mass difference is still negligible (between 0.003 and 0.0005 %).For the higher-resolution models there is also no systematic mass variation visible.We therefore can state that our implementation of the level set method is mass conservative.

Post-glacial rebound
Crameri et al. ( 2012) concluded that for a correct modeling of Earth's free surface by means of the sticky air approach the sticky air layer should either have a 5 orders of magnitude viscosity contrast with the underlying material and be 100 km thick or have a 4 orders viscosity contrast and be 200 km thick.In Fig. 4a we show results of a model which has a 100 km thick layer and 5 orders of magnitude difference in viscosity (η air = 10 18 Pa s) with the lithosphere.The gray area illustrates that the topography error of our model results never exceed 100 m.These results illustrate to our knowledge the first combination of the level set method and the sticky air approach.This result does require high local mesh resolution on the order of 1 km or less near the air-lithosphere interface.In order to minimize CPU time, and yet maintain a small resolution around the lithosphere-sticky-air interface, we investigated the possibility of reducing the thickness of the sticky air layer.This builds on the assumption that a certain thickness of the sticky air layer is required to ensure that the return flow in the air does not exert stress on the sticky air-lithosphere interface.The return flow is the result of the free slip top boundary condition (illustrated in Fig. 5a).Opening this top boundary, so making it a zero stress boundary which allows for through flow, makes the return flow within the sticky air layer disappear and strongly limiting the stress exerted on the interface (Fig. 5b).In Fig. 4b  and c

Subduction benchmark
In this benchmark the partly subducted slab retreats and then sinks into the mantle due to its higher density than the mantle.There is no overriding plate.For comparison we look at the slab tip depth over time and compare the results of our model with those of six model runs from Schmeling et al. (2008) (Fig. 6a).These six models are the fastest and slowest models of three different viscosity averaging schemes for tracers (harmonic, geometric, arithmetic) discussed in detail in Schmeling et al. (2008).Our diffuse boundary method (see Sect. 2) for viscosity can best be compared with the geometric averaging method for tracers.Figure 6a shows that our slab sinks a little slower than the geometric averaging models.The averaging of viscosity is especially important for large viscosity contrasts of which there are two in the models of Schmeling et al. (2008): the air-slab interface and the small entrained air layer in the subduction zone (Figs.7b and 9 of Schmeling et al., 2008).Both zones are important for the subduction velocity, the air-slab interface determines the decoupling between the two zones while a small entrained layer of air has a lubricating effect in the subduction zone.In our model the decoupling between slab and air is the same, but Fig. 7a illustrates that in our model due to the use of the level set method there is no entrainment of air and therefore no artificial lubrication of the subduction zone.This explains why the sinking slab in our models is slightly slower than the geometric averaging model results from Schmeling et al. (2008).

Slab detachment
The target of the slab detachment benchmark is the timing and depth of slab detachment through viscous necking.
Our results from two model runs with different mesh resolution are compared with results from the "V ∼ 100, top layer model" from Schmalholz (2011).There is ∼ 2 orders of magnitude difference in viscosity between the mantle and the location where necking occurs at startup. Figure 8 demonstrates resolution independence of our results and good agreement with the results from Schmalholz (2011).A particular benefit of using the level set method in monitoring  slab necking is that the moment the zero level set splits into two, disjoint domains can be used as a determination of the time of detachment.This was not determined by Schmalholz (2011).Figure 9 shows the time evolution of the necking process where acceleration of the process can readily be observed.In the first 17 Myr roughly half of the necking occurs while in the next 5 Myr it is completely detached.In the figures the thick white line represents the zero isocontour of the level set function.It is important to note that although the lithosphere has broken into two disjoint domains, it is still described by a single level set function.Using the data from Fig. 8 we can time the moment of detachment at 20.4 Myr.

Discussion and conclusions
The level set method is not often used in geodynamical modeling.To investigate its applicability and benefits we have conducted four benchmarks highlighting different aspects of geodynamical multi-material flow modeling.
All the benchmarks show the accuracy of the level set method.Our results for the Rayleigh-Taylor instability benchmark agree well with results published by other groups using other methods.When compared to earlier level set method results from Bourgouin et al. (2006) and Samuel and Evonuk (2010) our results agree better with the original paper of Van Keken et al. (1997).With this benchmark we also demonstrate that the level set method is mass conservative.The accuracy of our method can also be observed in the subduction benchmark and the slab detachment benchmark as our results in these cases also match previously published results.For the post-glacial rebound benchmark our results show an error of maximum 100 m with respect to the semi-analytical solution.This means we can resolve the long-scale topography well within 10 % of our local finiteelement resolution.All benchmarks thus demonstrate the level set method to be accurate.
In several cases we demonstrate favorable properties of the level set method compared to tracer-based methods.In the subduction benchmark we have shown that the level set method prevents entrainment of air into the subduction channel which otherwise would artificially lubricate the subduction zone.In the slab detachment benchmark we have demonstrated that the level set method can split from one bounded region into two bounded regions and accurately record the moment of detachment.In the post-glacial rebound benchmark we illustrated that we match the semianalytical solution using similar elemental resolutions as tracer-based methods.This further illustrates the statement of Zlotnik et al. (2008) that the level set method is favorable, particularly for 3-D applications, to the tracer method with regard to computational costs.
The level set function is chosen to be a signed distance function.This makes the implementation of a diffuse boundary method (Eq.8) simple and the boundary width (h) easily adjustable.In the Rayleigh-Taylor instability benchmark we have shown that such a diffuse boundary method has a stabilizing effect on the results but does not alter them in any significant way.
With the post-glacial rebound benchmark we demonstrated that the thickness of the sticky air layer can be reduced significantly when using a zero stress, free throughflow open top boundary resulting in an even better fit with the semi-analytical solution.
Overall we have shown that the level set method performs well and occasionally even better in geodynamical multimaterial flow benchmarks and could therefore be considered as an alternative for tracer-based and phase field methods.For 3-D applications one can add to this the lower computational costs compared to tracer-based methods (Zlotnik et al., 2008).
et al.: Geodynamical modeling using the level set method

Figure 1 .
Figure 1.(a) Setup for the Rayleigh-Taylor instability benchmark: material 1 has η = 100 Pa s and ρ = 1010 kg m −3 and material 2 has η = 100 Pa s and ρ = 1000 kg m −3 .(b) Setup for the subduction benchmark: material 1 is a sticky air layer with η = 10 18 Pa s and ρ = 0 kg m −3 , material 2 is a slab with η = 10 23 Pa s and ρ = 3300 kg m −3 and material 3 is the mantle with η = 10 21 Pa s and ρ = 3200 kg m −3 .(c) Setup for the post-glacial rebound benchmark: material 1 is a sticky air layer with η = 10 18 Pa s and ρ = 0 kg m −3 , material 2 is a lithosphere with η = 10 23 Pa s and ρ = 3300 kg m −3 and material 3 is the mantle with η = 10 22 Pa s and ρ = 3300 kg m −3 .(d) Setup for the detachment benchmark.Material 1 is a slab with nonlinear rheology and ρ = 3300 kg m −3 and material 2 is the mantle with η = 10 21 Pa s and ρ = 3150 kg m −3 .

Figure 2 .
Figure 2. Snapshots of the evolution of the density field and the level set field with time.Level set isocontours are plotted every 0.1.The thick white line represents 0.
Figure 3. (a) Root mean square velocity of the entire domain.The Van Keken data is the data from the marker chain method of Van Keken from the original article.(b) The relative mass difference over time according to Eq. (12).(c) A zoom-in of the first peak in the rms velocity plot of (a).

Figure 5 .
Figure 4. (a) The topography over time of our best fitting model run.(b) The data of open vs. free slip top boundaries for different sized sticky air layers.The solid green, purple and blue line exactly overlap.Confirmed by the close-up in (c).(c) Contains a close-up of (b) showing that the three open boundary model runs (solid lines) do indeed yield exactly the same result.

Figure 6 .Figure 7 .
Figure 6.(a) Depth of slab tip vs. time.(b) through (e) snapshots in time of the slab.The red line indicates the location of the zero isocontour of the level set function tracking the surface and the black line indicates the zero isocontour of the level set function tracking the slab.

Table 1 .
Table containing the timing and hight (max v rms ) of the first peak of our model runs and selected others from the literature.