New analogue materials for nonlinear lithosphere rheology, with an application to slab break-off

Abstract Stress-dependent nonlinear upper mantle rheology has a firm base in rock mechanical tests, where this nonlinearity results from dislocation creep of minerals. In the last few decades there has been some attention to nonlinear, power-law, materials for application in scaled analogue experiments for tectonic processes. However, studies describing the rheology of analogue materials with the same nonlinear dependency on stress as observed for lithospheric mantle materials at relevant stress levels, are still lacking. In this study we have developed and rheologically tested materials based on combinations of silicone polymers and plasticine, with the aim of obtaining a material that can serve as a laboratory analogue to the power-law rheology of olivine aggregates at lithospheric mantle conditions. From our steady-state creep tests we find that it is possible to obtain such a power-law material, with effective viscosities over relevant model stress ranges [5–4000 Pa] that allow for nonlinear deformation at laboratory time scales. We apply the developed material to a process where localized deformation of the lithosphere can be expected: slab break-off. We study this process using analogue models, where we apply the new nonlinear material to the lithospheric mantle domains, while we use Newtonian glucose to represent the low viscous asthenosphere. Now that we properly manage power-law behavior in our analogue lithosphere materials, we are able to model localized lithospheric tearing.

In this study we have developed and rheologically tested materials based on combinations of silicone polymers and plasticine, with the aim of obtaining a material that can serve as a laboratory analogue to the power-law rheology of olivine aggregates at lithospheric mantle conditions. From our steadystate creep tests we find that it is possible to obtain such a power-law material, with effective viscosities over relevant model stress ranges  * Email: d.b.t.broerse@uu.nl; Corresponding author that allow for nonlinear deformation at laboratory time scales.
We apply the developed material to a process where localized deformation of the lithosphere can be expected: slab break-off. We study this process using analogue models, where we apply the new nonlinear material to the lithospheric mantle domains, while we use Newtonian glucose to represent the low viscous asthenosphere. Now that we properly manage power-law behavior in our analogue lithosphere materials, we are able to model localized lithospheric tearing.

Introduction
Permanent deformation of the Earth's lithosphere can occur in brittle or ductile (viscous) modes, where the latter prevails for increasing temperatures and depths (Goetze and Evans, 1979;Kohlstedt et al., 1995). Ductile deformation involves a permanent change of shape without fracturing, and can be caused by various microphysical creep mech-anisms, each of which results in deformation rates with different dependencies on external conditions such as stress or temperature, and microstructure of the material (Ranalli, 1995). For ductile steady-state deformation we can discriminate between Newtonian rheology, where strain rate is proportional to stress and hence viscosity is constant, and nonlinear rheologies, where viscosity has a nonlinear dependency on stress (Ranalli, 1995). Nonlinear ductile deformation often comes in the form of powerlaw creep, where strain rate˙ has a power-law relation to stress σ, leading to material weakening at elevated stresses, a phenomenon often called pseudoplasticity:˙ = Cσ n n > 1 (1) Many different materials, from tomato puree to metals exhibit both types of ductile deformation: linear Newtonian creep at low stresses and powerlaw creep at high stresses (Barnes, 1999). Also for mantle materials, such as olivine aggregates, laboratory tests at high temperatures show this transition from Newtonian to power-law creep with increasing differential stress (Hirth and Kohlstedt, 2003).
The existence of plate boundaries implies localized deformation, not only in the brittle, but also ductile regime in relatively weak shear zones (Bercovici, 1993). Localization of ductile deformation is often attributed to strain softening as a result of previously accumulated deformation (e.g. Poirier (1980); Bercovici and Karato (2002)), by mechanisms such as recrystallization (Karato et al., 1980) or shear heating (Regenauer-Lieb and Yuen, 1998). The amount of weakening or strengthening due to strain softening or hardening, respectively, in shear zones is however debated (De Bresser et al., 2001;Drury, 2005;Platt and Behr, 2011;Tasaka et al., 2017). Alternatively, in settings with highly varying stresses, power-law creep may result in high rates of deformation, as the strain rate will localize stronger than the stress itself. Such a system weakening will be especially efficient when it results in geometrical changes of tectonic units such that stress becomes increasingly localized by a positive feedback, as is observed in slab necking or tearing processes in numerical and analytical models (Andrews and Billen, 2009;Schmalholz, 2011). Several numerical studies suggest the importance of power-law stress weakening of the lithosphere for transient localization phenomena; in subduction models power-law creep has been invoked as mechanism required for slab break-off (Davies and von Blanckenburg, 1995;Van Hunen and Allen, 2011;Duretz et al., 2011). Furthermore, power-law creep has been invoked to explain rapidly decaying postseismic deformation following several earthquakes (Freed et al., 2006(Freed et al., , 2012Masuti et al., 2016;Sobolev and Muldashev, 2017). In combination with strain softening a strongly stress-dependent rheology is a precondition for stable shear zones (Rutter, 1999).
Analogue models of tectonic processes aim at constructing a scaled representation of a natural tectonic setting, where application of geometrical, kinematic and dynamic scaling of the natural setting to the laboratory analogue is required (Hubbert, 1937;Weijermars and Schmeling, 1986). A precondition for correct dynamic similarity is that the rheological properties of the materials used are properly scaled with respect to the rheological properties of the respective natural rocks. Often in analogue models of tectonic processes, only the Newtonian behavior of rocks deforming in the ductile regime is modeled, by using materials that show near-Newtonian behavior, such as silicone polymers (Weijermars, 1986;Ten Grotenhuis et al., 2002;Rudolf et al., 2016). While for many tectonic settings the assumption of Newtonian behavior can be justified, processes that may involve a high variability of differential stresses call for materials that exhibit both Newtonian and power-law rheology.
Over the last few decades there have been a number of studies that tested power-law materials for the application in analogue models of tectonic processes (see an overview in section 4). However, rheological studies of analogue materials that are tailored to represent the rheology of the lithospheric mantle are currently missing. Especially, there is a lack of materials that show the same nonlinearity to stress, expressed by the stress exponent n, as the most abundant upper mantle material olivine (n ≈ 3.5 Hirth and Kohlstedt (2015)). For existing analogue materials, the power-law nonlinearity is either too weak, or overly strong compared to olivine aggregates (see section 4.2). Furthermore, rheological data for power-law lithosphere analogue materials is lacking at the low stresses that are typical for analogue models at the upper mantle scale (< 100 Pa, Boutelier and Oncken (2011)). In our paper we develop materials based on silicones and plasticine mixtures -informed by previous works of Weijermars (1986); Sokoutis (1987); Zulauf and Zulauf (2004) and Boutelier et al. (2008) -and study their behavior to determine how these can serve as mechanical analogs for nonlinear ductile creep in the lithospheric mantle. Furthermore, we are interested in the type of deformation behavior of the developed materials, distributed or localized, in tectonic settings where lithospheric tearing may occur.
This paper consists of two parts: in the first part we estimate the average rheology of lithospheric mantle and test the rheology of silicone-plasticine materials at stress levels typical for analogue models at the mantle-lithosphere scale. In the second part we present an application in a model of break-off of subducted lithospheric slabs, where localisation of strain is a prerequisite and a powerlaw material may potentially facilitate this mode of deformation.
2 Nonlinear rheology of the lithospheric mantle

Laboratory flow laws
We intent to design a material that can function as a rheological analogue for the average lithospheric mantle rheology. Olivine ((M g, F e) 2 SiO 4 ) makes up roughly half of the upper mantle composition at depths down to 400 km (Putnis, 1992) and is the weakest upper mantle constituent (Karato and Wu, 1993). Therefore, olivine aggregate rheology is generally thought to represent the rheology of the upper mantle (e.g. Goetze and Evans (1979); Hirth and Kohlstedt (2003)). The steady-state rheology of olivine aggregates at high temperatures shows a transition from a viscosity that is stressindependent to a stress-dependent viscosity at increasing stress (Mei and Kohlstedt, 2000b). After application of a constant stress, olivine shows recoverable deformation, both immediate (elastic) and time-dependent (anelastic), next to permanent (viscous) deformation. After a transitional period where both anelastic and viscous deformation occur, strain rates decrease and converge to a steady-state regime (Chopra, 1997). Steady-state strain rates from rheological experiments on olivine can be explained from contributions from diffusion creep (Newtonian) and dislocation creep (powerlaw) (Hirth and Kohlstedt, 2003): where power-exponent n has been estimated to be in the range 3 -4.9 (Mei and Kohlstedt, 2000a;Korenaga and Karato, 2008) and is likely close to 3.5 (Hirth and Kohlstedt, 2015). Parameters C are functions of grain size, temperature, pressure and water content (Hirth and Kohlstedt, 2003): with material parameter A, grain size d, grain size exponent p, C OH the water content with exponent r, activation energy E, activation volume V , pressure P , gas constant R and temperature T in Kelvin.

Ellis rheology
Comparison and scaling of rheologies is simpler when the strain-rate vs. stress equation (2) is rewritten to a viscosity-stress relation: (4) which we rewrite to the Ellis rheological model: with: The Ellis model describes the combination of power-law and Newtonian rheological behavior in terms of the zero-stress Newtonian viscosity η 0 and transition stress σ t (Barnes, 1999). At the transition stress, the strain rate contributions from the Newtonian and nonlinear mechanisms in equation (2) are equal, and σ t thus marks the change of rheological behavior: Newtonian for stresses below σ t and power-law above σ t .

Average ductile lithosphere rheology
We are interested in an average ductile lithospheric mantle rheology, which we cast in terms of an Ellis rheology. For first order estimates of the variability of η 0 and σ t , we use equations (3),(6) and (7), where we apply parameters for dry and wet (constant C OH ) olivine diffusion and dislocation creep from Hirth and Kohlstedt (2003), table 1, including the value for n = 3.5 for dislocation creep. We apply a correction to C for parameters based on uni-axial experiments: C = C 2 3 (n+1)/2 (Ranalli, 1995). We assume a linear geotherm between the Moho and 1350 • C at the base of the thermal boundary layer, where for continental lithosphere we place the Moho with a 400 • C temperature at 40 km depth and the 1350 • C isotherm at 200 km (e.g. Jaupart and Mareschal (2007)); and for an oceanic lithosphere we situate the Moho with a 50 • C temperature at 5 km depth, and the 1350 • C isotherm at 100 km depth. We define the ductile lithosphere between a brittle-plastic transition at a strain rate of 10 −14 s −1 using the Goetze criterion for the brittle strength of the lithosphere (σ 1 − σ 3 < P ) (Kohlstedt et al., 1995) and a lithospheric-asthenospheric boundary (LAB) defined by η 0 = 10 20 Pa s.
Using the described wet and dry olivine rheologies, the continental and oceanic geotherms and a grain size range between 0.1 and 10 mm, we estimate η 0 and σ t as a function of depth. The first two panels of figure 1 show results for an oceanic lithosphere with a dry rheology, where we normalise η 0 by η 0,LAB , the diffusion creep viscosity that defines our LAB. Figure 23 (B) shows similar results for η 0 and σ t for a wet olivine rheology, respectively continental lithosphere. The zero-stress viscosity η 0 of the lithospheric mantle, for a constant grain size throughout the lithosphere, is highly variable, mainly due to the strong temperature dependency of diffusion creep. For lithospheric temperatures, the transition stress σ t is dominantly dependent on the grain size due to the grain size dependency of diffusion creep: σ t increases for decreasing grain size. For the last panel of figure 1 we take the average of η 0 , and use σ t from the uppermost part of the ductile lithosphere (as the most viscous part should be representative for the overall strength) to arrive at an average rheology. We repeat the same analysis for wet as well as dry rheologies and continental as well as oceanic lithospheres. While for a given pressure and temperature the olivine creep law consistently results in a lower viscosity for a wet rheology compared to a dry rheology (figure 23), this lower viscosity also results in a shallower brittle-plastic transition, resulting in a somewhat higher average η 0 . Figure 1 shows that the average η 0 is not well constrained, as values are very sensitive the choice of grain size and to a lesser extent to the presence of water and type of lithosphere. Average ratio η 0 /η 0,LAB ranges between 2 · 10 2 and 4 · 10 8 , with the lower value corresponding to the smallest considered grain size (0.1 mm). Similar large increases in viscosity throughout the lithosphere can be found in figure 2 of Karato and Wu (1993) Figure 1: Average lithosphere rheology using olivine rheological parameters from Hirth and Kohlstedt (2003), left panel: η0 for an oceanic lithosphere and a dry rheology for different values of grain size d, dotted lines denote average η0; middle panel: σt for an oceanic lithosphere and a dry rheology, dotted lines denote values below the lithospheric-asthenospheric boundary (LAB) defined by η0 = 10 20 Pa s. The top depth represents the brittle-plastic transition defined by the intersection of the Goetze criterion and the ductile stress at a strain rate of 10 −14 s −1 . Right panel: combined results for oceanic and continental lithospheres with wet or dry olivine rheologies and 0.1 and 10 mm grain sizes, using the average η0 and maximum σt and n = 3.5. or figure 2.1 of Freeburn (2016); Billen and Hirth (2007) consider a 10 6 − 10 7 viscosity contrast for a subducting slab with respect to the asthenosphere, based on the same rheological parameters by Hirth and Kohlstedt (2003). The values for σ t span a wide range: 2 · 10 6 − 3 · 10 9 Pa, comparable to the transition stresses in figure 6 from Drury (2005). Nevertheless, once the power-law behavior is dominant at higher stresses, the right panel of figure  1 shows that the rheology follows a narrow band for all tested parameters.
Regarding parameters that describe recoverable deformation, the elastic shear modulus µ e for the mantle lithosphere is on the order of 67 GPa (Dziewonski and Anderson, 1981), and anelastic olivine shear moduli µ a have been estimated to be in the range 2.3 -37 GPa (Chopra, 1997). We have thus obtained ranges for the three parameters of the Ellis model (n, σ t and η 0 ) that can describe the average lithospheric mantle rheology, next to shear moduli that constrain the amount of recoverable deforma-tion. We will use these values in the next sections for developing an analogue lithospheric material.

Dynamic scaling of non-Newtonian materials
Weijermars and Schmeling (1986) developed rules for dynamic similarity of power-law materials in the case of negligible inertia. We propose a simpler approach using the Ellis rheological model (equation (5)) for experiments performed under standard gravity. The Ellis model equation contains three variables for which we can make use of existing scaling laws. Namely, the stress exponent n should be equal for model analogue (m) and natural (n) materials (Weijermars and Schmeling, 1986): Furthermore, similarity of viscous forces requires constant viscosity ratios between model and nature values (Hubbert, 1937), which implies that the zero-stress viscosity η 0 scales with the same ratio as all other viscosities: Finally, for σ t we can simply use the existing scaling relations for stresses (e.g. Weijermars and Schmeling (1986)): where ρ is density and l is length. Alternatively, if dynamic scaling is performed by using density differences between the sinking material and the surrounding fluid ∆ρ (e.g. Ribe and Davaille (2013); Schellart and Strak (2016)), the expression becomes: The expression for transition stress σ t implies that model dimension scale is important in the transition from linear to nonlinear behavior for a given material rheology. Time furthermore scales as (Ramberg, 1981;Schellart and Strak, 2016): 4 Nonlinear analogue materials

Required analogue material rheology
Analogous to olivine aggregate rheology (section 2), the required lithospheric mantle material should behave Newtonian at low stresses, with a transition to power-law creep at increasing stress, with a stress exponent around 3.5. This parameter and other relevant rheological parameters for the analogue material are summarized in table 1. The table shows that the zero-stress viscosity η 0 for the lithospheric mantle should be in the range 2 · 10 4 -4 · 10 10 Pa s (when the asthenospheric mantle is modeled using glucose), and that the transition stress σ t should fall somewhere between 0.2-300 Pa, depending on the size of the model and the scaling of densities, as explained by the caption of table 1. Values for η 0 and σ t are highly correlated: large σ t apply to low values of η 0 . Elastic and anelastic strains in settings with large permanent deformation should be small, a constraint which is satisfied when the shear moduli (approximately) scale to the natural values. The shear modulus is scaled in the same manner as stress (Hubbert, 1937) (using equation (11) or (10)), which implies that for the analogue material the elastic shear modulus should be on the order of 10 4 Pa, and the anelastic shear modulus in the range 10 3 -10 4 Pa.

Previous studies on power-law analogue materials
Over the last few decades there have been a number of studies that apply power-law (pseudoplastic) materials in analogue models of tectonic processes. Figure 2 summarizes the rheology of existing materials with stress dependent ductile rheologies. We also include commonly used Newtonian silicone polymers for reference such as PDMS (Rudolf et al., 2016) and rhodorsil gum (Treagus and Sokoutis, 1992). The required power-law behavior is indicated in figure 2 by the diagonal lines that represent a stress-exponent n of 3.5. The grey curved band represent the appropriate viscosities for analogue lithospheric mantle materials as a function of applied stress, using the requirements from section 4.1 (figure 1) in combination with an analogue asthenospheric viscosity of ≈ 10 2 Pa s (i.e., glucose). Figure 2 shows that none of the available materials can be used as a lithospheric mantle analogue without any modifications or additional rheometric measurements. Plasticines, either pure (McClay, 1976;Treagus and Sokoutis, 1992;Barnes, 1999) or mixed with an oil (Zulauf and Zulauf, 2004) have too high viscosities, while paraffin has a similar disadvantage (Rossetti et al., 1999). Furthermore plasticines weaken too drastically for elevated stress, indicated by slopes that are steeper than the diagonals that represent n = 3.5; for plasticine stress exponents n have been estimated to be in the range 5-10 at room temperature (Weijermars, variable natural values scale factor n m target analogue values n 3.5 1 3.5 η 0,lith /η 0,asth 2 · 10 2 − 4 · 10 8 1 2 · 10 2 − 4 · 10 8 η 0,lith [Pa s] 2 · 10 22 − 4 · 10 28 ≈ 10 18 2 · 10 4 − 4 · 10 10 σt [Pa] 2 · 10 6 − 3 · 10 9 10 6(I) − 10 7(II) 2 − 300 (I) /0.2 − 30 (II) µe [Pa] 7 · 10 10 10 6(I) − 10 7(II) 7 · 10 4(I)/3(II) µa [Pa] [0.2 − 4] · 10 10 10 6(I) − 10 7(II) [0.2 − 4] · 10 4(I)/3(II) prescribes that stress exponent n should be equal for nature and model. The scaling of stresses and moduli is dependent on the length scale of the model, as well as density or density difference scaling, which means that there is a possible range for the scaling of stress. We present two end-member scaling options; scaling option I uses a density difference scaling of ∆ρn ∆ρm = 0.5 and uses a model height of 20 cm to represent an upper mantle depth of 660 km. This gives a stress scale σn/σm ≈ 10 6 (equation (11)). This option resembles the setup of our models in the second part of this paper. For scaling option II we assume a 10 cm high model and a density ratio ρn ρm = 2 (in line with densities of commonly used analogue materials), resulting in σn/σm ≈ 10 7 (equation (10). From equation 9 follows that the ratioη 0,lith /η 0,asth should be equal for nature and model. The value for the average analogue lithospheric viscosityη 0,lith depends on the model asthenospheric viscosity, in our case 2 · 10 4 <η 0,lith < 4 · 10 10 Pa s when η asth ≈ 10 2 Pa s (i.e., glucose, Schellart (2011)).
1986; Schöpfer and Zulauf, 2002;Zulauf and Zulauf, 2004). However, an advantage of plasticine is its Newtonian behavior at low stresses, roughly below a transition stress σ t of 10 4 Pa (Barnes, 1999). The rheology of the silicone putty Dow Corning 3179 falls within the indicated area, but the nonlinearity of the material is too high, expressed by a stress exponent of ≈ 6 (Treagus and Sokoutis, 1992). The petrolatum + paraffin oil mixtures (designed to represent a subduction interface (Duarte et al., 2013)) display a reasonable Newtonian platform at low stresses as well, however these materials suffer from the same extreme weakening in the nonlinear regime as plasticines, which seems approximately constant for different petrolatum to paraffin oil ratios.
Carbopol has overall lower viscosities, but has been described to exhibit a yield stress (i.e., a stress below which no flow occurs) instead of a low-stress Newtonian plateau (Di Giuseppe et al., 2015), and viscosity decreases excessively at the stress range of interest. Natrosol solutions such as described by  nicely show a constant viscosity for low stresses, a slightly too low stress exponent (n ≈ 2), while the viscosities are less than needed for lithospheric application.
The silicone polymer (PDMS) -plasticine mixtures as tested by Boutelier et al. (2008) approach the target rheology, but have either a slightly too low stress exponent (n = 2.8) for 20 volume % of plasticine or too high overall viscosity (n = 4.6) for 40 volume % of plasticine. Next, the stress levels over which the rheological parameters have been obtained are rather limited, and do not cover the expected stress range for upper mantle laboratory models. Silicone polymer -plasticine mixtures seem the most promising materials, which however need further adjustments to conform to the target rheology over wider stress ranges.

Material use in this study
Our goal is to develop an analogous material for olivine aggregate rheology at lithospheric conditions, combining both linear viscous (diffusion creep) behavior as well as power-law behavior, with a stress exponent n in the range 3-4, a zerostress viscosity η 0 in the range 2 · 10 4 -4 · 10 10 Pa s and a transition stress σ t in the range 2-300 Pa (corresponding to the light grey area in figure 2). We choose to work with plasticine-PDMS mixtures, as Boutelier et al. (2008) show that the stress exponent can be adjusted by the ratio of PDMS to plasticine. Plasticine-silicone polymer materials have been used or tested in previous studies (Weijermars, 1986;Sokoutis, 1987;Treagus and Sokoutis, 1992;Schrank et al., 2008), and we aim to overcome the problem of the high viscosity of these materials by adding a low viscosity silicone oil, which is likely to have only limited effects on the stress exponent but can significantly reduce the overall viscosity (Zulauf and Zulauf, 2004). Both plasticine and PDMS exhibit Newtonian behavior at low stresses in accordance with the Ellis rheology. We use the following material components to develop a lithospheric mantle analogue material:

Plasticine
Plasticine, and plasticine-like materials made by various manufacturers, are non-drying materials, and sold commercially as modeling clays. The components of plasticine or its varieties that go by names such as plastilin or plastilina are usually confidential, but generally these materials consist of an organic matrix (wax, vaseline, oils) with inorganic fillers (McClay, 1976). Creep tests have shown the limited applicability of many plasticines, which often do not reach steady-state creep but exhibit increasing viscosities with increasing strain (Zulauf and Zulauf, 2004;Boutelier et al., 2008). Zulauf and Zulauf (2004) describe a plasticine-like material, that they label Beck's orange, containing organic fillers (potato starch) and that shows only limited strain dependency. However, this material is no longer in production. For our experiments we use a similar material that is currently produced by the manufacturer Becks Plastilin and sold under the name Natur-Kreativknete, and that we refer to as organic plasticine. Relevant properties are listed in table 2. According to the manufacturer the organic plasticine is composed of beeswax, natural oils and plant-based filler (potato starch). We observe that plasticine is not very self-adhesive, which means that seams originating from the assemblage of different tectonic units may act as weak zones.

PDMS
The silicone polymer polydimethylsiloxane, PDMS, is a commonly used material in analogue models, especially for its Newtonian rheological properties for a wide range of strain rates (Johnson, 1961;Weijermars, 1986). PDMS is produced in a continuous range of viscosities, where longer polymer lengths result in higher viscosities (Weijermars, 1986). The high-viscosity PDMS that we use in our study, SGM36 produced by Dow Corning, has a viscosity of about 2·10 4 Pa s for stresses below 4 kPa, or alternatively, below strain rates below 2 · 10 −1 s −1 (Rudolf et al., 2016). For higher stresses, or strain rates, the material increasingly starts to behave like a power-law material. One of the advantages of PDMS is the self-adhesive character, allowing the construction of multi-layered tectonic models without weak zones stemming from the assembly of the model.

Low-viscosity silicone oil
Zulauf and Zulauf (2004) have shown that addition of oil to plasticines reduces the effective viscosity while affecting the stress exponent only in a minor way. We add various amounts of silicone oil (50 cs Xiameter PMX-200 fluid), which is chemically similar to the SGM36 PDMS, but due its shorter polymer chains has a viscosity six orders of magnitude lower than PDMS (see table 2).

Iron powder
PDMS, low viscosity silicone oil and organic plasticine have similar densities, all about 450 kg m −3 lower than the density of glucose (table 2), the material that we use to represent the low-viscosity sub-lithospheric mantle (e.g. Davy and Cobbold (1988); Faccenna et al. (1999); Guillaume et al. (2009);Duarte et al. (2013)). By adding fine iron powder (Funiciello et al., 2004;Schellart, 2008; sample h R σ r Newtonian Power-law rs Figure 3: Schematic rheometer setup, with gap height h and disk radius R, including a schematic depiction of shear stress σ variation as function of radial distance r for Newtonian and power-law materials, rs is the radial distance where the shear stress curves of Newtonian as well as power-law materials intersect (Carvalho et al., 1994). Guillaume et al., 2013) we adjust the mixture density to be able to create both buoyant and negatively buoyant materials. We mix and homogenise materials by hand. Table 3 lists the composite materials tested in this study.

Measurement setup
We use an AR-G2 rheometer (by TA Instruments) that we operate in a parallel plate setup, as shown in figure 3. In all measurements the radius R of the upper plate is 12.425 mm and a gap height of 2.5 mm is applied. The rheometer applies a torque M and measures the angular deflection θ. Strain γ (engineering strain γ = 2 ) is a simple linear function of radial distance r: Whereas the applied torque M is the integral of the product of shear stress and arm over the area: As Macosko (1994) shows, the shear stress at the rim (r = R) can be written as a function of the torque and the derivative of torque to strain: The derivative is unknown, except for Newtonian materials, where it is 1. For that reason the rheometer reports apparent shear stress σ a values at the rim that are equal to the true shear stress σ if the material would behave as a Newtonian material: To calculate the actual stress (for a non-Newtonian material) either the derivative of the torque to strain rate should be measured, or one can make use of the fact that at one radial distance (r s ) the curves for apparent Newtonian σ a and true power-law stresses σ intersect (see figure 3): As shown by Carvalho et al. (1994) at radial distance r s = 0.755 * R, the apparent and true shear stresses are equal with an error less than 1% for power-law materials with stress exponents n larger than 0.8. As both strain and apparent shear stress vary linearly with radial distance, we can apply this so-called one point correction by taking shear stress and strain at r s :   We test the materials under constant stress, with apparent stress levels σ a in the range [5;4000 Pa], in a creep and recovery setup. Because of our application to modeling tectonic deformation, we are interested in large and permanent deformation. A creep and recovery setup allows for application of large strains, and due to the relaxation phase after each loading phase, it can discriminate between recoverable strain and permanent strain, where it does not need a priori assumptions on the type of rheology for the calculation of strain rates. Before the first application of stress we let the sample rest in the rheometer for at least 3 hours, to relax residual stresses. Stresses are applied in ascending order; figure 4 shows a typical relaxation phase after a period of loading. The duration of each stress application and subsequent relaxation differs for each stress level, where we aim at creating permanent deformation during stress application and (nearly) full relaxation during the rest period.

Data analysis
The creep and recovery tests show that the recorded strain contains a considerable elastic and anelastic component; from figure 4 we can observe immediate elastic (recoverable) deformation, followed by a superposition of a transient anelastic (recoverable) and steady-state (permanent) deformation. Especially for low stresses, where long periods of loading (several hours) are needed to produce significant permanent deformation, it is difficult to discriminate between the early transient deformation and the permanent deformation, which deviates slightly from a constant strain rate due to strain hardening or strain softening. While the strain curve resembles the theoretical curve for a Burgers body at constant stress (see Findley et al. (1976)), observed strain hardening (at low stress) or strain softening (at high stress) does not permit the application of the Burgers model.
In the absence of a mathematical model that would allow us to estimate a steady-state viscosity at each stress level, we therefore determine the mean steady-state viscosity per stress level. The total strain that has accumulated during the loading period up to the unloading time t u comprises of elastic γ e , anelastic γ a and permanent γ p strains, where γ e and γ a recover during the relaxation phase, such that the remaining strain equals to γ p (figure 4). The average steady-state viscosity then becomes: To describe the stress-dependent viscosity we fit the Ellis model of equation (5) to the average viscosityη as a function of stress σ, or solely a powerlaw equation (equation (1)) in case no transition from linear to nonlinear viscosity could be observed. We perform the fitting in a weighted nonlinear least squares sense using the Levenberg-Marquardt algorithm, where we minimize: using weights W = 1/η(σ).
For strain rates below 10 −4 1/s, the rheometer reports increasingly noisy strain rates. To determine instantaneous effective viscosities η eff , we therefore first smooth the strains using smoothing splines, and afterwards differentiate the smoothed strainsγ with respect to time to create smoothed strain rateŝ At each stress level we derive the elastic µ e and anelastic shear modulus µ a from the elastic and anelastic strains (see figure 4) determined from the relaxation curves: Effective viscosity shows a relatively large increase during the start of the loading, which is most pronounced at lower stresses. In this phase the strain rates are dominated by transient, anelastic deformation. For those stresses where the average steady-state viscosityη is close to the observed effective viscosityη eff (depicted by a dot) the transient phase is over. We can still observe small (on logarithmic scale) variations in viscosity in the phase when only permanent deformation is occurring; here strain softening (decreasing viscosity) is observed for high stresses, strain hardening (increasing viscosity) for low stresses. For curves shown in grey there was an unknown amount of permanent deformation, as during the end of the chosen period for relaxation of the sample the strain rate was still too high (more than 15% of the final loading strain rate). All loading and relaxation strain curves are included in appendix A, whereas the rheometry data are published in a separate data publication (Broerse et al., 2018).
As shown by figure 6, five of the six materials (all except the sample 80OPL 20PDMS) exhibit a clear nonlinear stress dependent viscosity in the applied stress ranges. The material with the lowest plasticine fraction (40OPL 60PDMS) shows the lowest zero-stress viscosity η 0 and the weakest stress nonlinearity, demonstrated by the power exponent n of 2. Increasing the plasticine content increases the power exponent, as expected: all materials having a plasticine to PDMS mass ratio of 60:40 have power exponents in the 3-4 range. Only for higher silicone oil contents (24% and 16% mass fractions) we are able to observe a transition from Newtonian to nonlinear steady state viscosity. The addition of the low-viscosity silicone oil lowers the viscosity over the whole stress range, but seems to have a minor influence on the power exponent n (comparable to the findings of Zulauf and Zulauf (2004) Figure 5: Effective smoothed viscosityη eff versus time for the six tested materials. The coloured dot denotes the average viscosity based on the observed permanent strainη. Grey lines stand for measurements for which we were not able to extract the permanent strain with sufficient confidence (i.e., the strain rate at the end of the relaxation period was higher than 15% of the strain rate at the end of the loading period). Numbers stand for stress level in Pa.
power-law equation to the curves of these materials). At even higher fractions of plasticine (80%) we obtain only a weak stress-dependence of the viscosity. As the parallel-plate setup leads to slip between material and plates at stresses at the high end of our applied stress range, we are not able to determine a possible transition to nonlinear viscosities for this material at higher stresses. Table 4 summarizes the estimated parameters, next to the confidence intervals.

Elastic parameters
As figure 7 shows, the elastic shear modulus γ e changes only little with stress, and decreases for larger fractions of PDMS and low viscosity silicone oil. Table 4 summarises the average shear modulus for each material. The anelastic shear modulus γ a is stress dependent: while constant at those stress intervals where the steady-state viscosity is Newtonian, γ a increases in the power-law viscosity stress domain, implying a decreasing amount of anelastic deformation at increasing stress.

Temperature effects
To infer the viscosity sensitivity to temperature we test the material 60OPL 40PDMS 45.6Fe 24oil at 30 • C. Figure 8 shows a non-uniform and modest reduction in viscosity (by a factor of two on average) for an increase in temperature from 20 • to 30 • C. The degree of nonlinearity, expressed by stressexponent n is not significantly affected, see table ??.

Iron filler effects
Up to this point we have adjusted the iron powder content to achieve a uniform density for all samples. For analogue models, variations in density are needed to accommodate naturally occurring density variations throughout the lithosphere. To that end we test the influence of significant iron filler content variations on the rheology. Figure 9 shows that large changes in iron powder contents, which we vary between 37 -92% in mass, have relatively small effects on the rheology. A higher iron filler content seems to slightly lower the viscosity at stresses below 100 Pa, and slightly increase viscosity above. The power-law behavior is hardly affected based on visual inspection of figure 9, but the estimated exponent n is lower for the 92.0 mass % iron sample. However, this can merely be attributed to a worse fit to the data, as the sharp kink in viscosity at around 20 Pa cannot be properly represented by the Ellis model.

Discussion of analogue material rheology
In our search for analogs for the ductile deformation of olivine aggregates we find that the materials with a plasticine to PDMS mass ratio of 60:40 have the desired rheological characteristics. Firstly, the material shows a nonlinear stress dependence (expressed by a power exponent of around 3.5) that nicely matches the stress dependence reported for olivine (in the range 3-4.9 e.g. (Mei and Kohlstedt, 2000a;Korenaga and Karato, 2008)). Secondly, the zero-stress viscosity of about 10 8 Pa s (for a 24% inclusion of silicone oil) allows a viscosity contrast with the asthenosphere in the intended range of 2 · 10 4 -4 · 10 10 Pa s when glucose is used to represent the asthenosphere (η a ≈ 10 2 Pa s). The material can thus represent the average viscosity of that part of the lithosphere that deforms in a ductile regime. Furthermore, the transition stress σ t (for 24% silicone oil added, see table 4) is in the appropriate range compared to scaled values for olivine: 2-300 Pa (see section 4). This indicates that the power-law behavior as witnessed in our materials can be effectuated in properly scaled laboratory models, as this transition occurs at realistic tectonic stresses. Finally, we observe that elastic and anelastic strains are limited and that both elastic as well as anelastic shear modulus properly scale with natural values (target γ e and γ a shear moduli are approximately 10 4 and 10 3 -10 4 Pa, respectively, section 4.1), which implies that the viscous behavior of the material is dominant for large deformation. All mixtures have self-adhesive  (table 3) with a fit (equation (6)) to the theoretical models: power-law, Ellis and Newtonian rheologies. The grey area represents the desired range in viscosity vs. stress as described in section 4 for stress scale option I. 2.8·10 4 4.6·10 3 -3.3·10 5 60OPL 40PDMS 45.6Fe 24oil 3.6(3.2 -4.1) 2·10 8 (4·10 7 -4·10 8 ) 18.5(7.6 -29.5) 2.4·10 −12 1.8·10 4 2.8·10 3 -5.2·10 5 80OPL 20PDMS 13.7Fe 1 2·10 9 (3·10 8 -3·10 9 ) --1.1·10 6 1.2·10 5 -4.4·10 5 Table 4: Estimated parameters of steady-state viscosity as a function of stress (see figure 6): stress exponent n, zero-stress viscosity η0, transition stress σt, power-law material constant C, average elastic modulus µe and range of anelastic shear modulus µa. Values within brackets provide the 95% confidence intervals.  properties, which prevents material boundaries to persist as weak zones in the final model assembly.
The occurrence of a transition from Newtonian behavior to power-law behavior, as described by Barnes (1999) has been subject of recent debate. The Newtonian plateau at low stress for Carbolpol solutions -a material that has also received attention to serve as an analogue model for tectonic studies, e.g. Schrank et al. (2008); Di Giuseppe et al.  (6)). The grey area represents the desired range in viscosity vs. stress as described in section 4 for stress scale option I. Figure 9: Steady state viscosities per stress level for 60OPL 40PDMS 24oil with three different iron filler ratio's, with an Ellis fit (equation (6)). The grey area represents the desired range in viscosity vs. stress as described in section 4 for stress scale option I.
(2015) and figure 2 -observed by Barnes (1999) has been interpreted by Møller et al. (2009) as an experimental artefact (see also Denn and Bonn (2011); Dinkgreve et al. (2017)). These studies show that at low stresses the effective viscosity of a Carbopol solution increases with time, which these authors attribute to anelastic deformation rather than Newtonian steady-state viscosity. They claim that there exists a yield stress below which no flow occurs, and that marks the transition material transforms from a solid to a liquid state. For two of our tested materials (40OPL 60PDMS and 60OPL 40PDMS 24oil) we find permanent strain for the complete applied stress range, which means that flow occurs down to the lowest stress level. For these two materials we observe a Newtonian plateau at low stress, a similar behavior as Boutelier et al. (2016) find for the low-viscosity natrosol solutions (η ≈ 10 1 − 10 4 Pa s, see figure 2). The Ellis model fits well to the observed combined low-stress Newtonian and high-stress power-law viscosities, only at the highest stress levels the fit becomes less and the observations suggest a decrease in stress exponent.

Material application: slab break-off
To demonstrate the capability of our developed power-law material to create localized deformation in settings with significant variation of differential stress and a positive stress feedback, we apply the developed material in a model of continental collision and possible slab break-off (Davies and von Blanckenburg, 1995;Wong A Ton and Wortel, 1997).
The subduction of oceanic lithosphere and the following closure of the ocean leads to the subsequent forced subduction of buoyant continental lithosphere. Whereas the negative buoyancy of oceanic lithosphere initially drives sinking of the oceanic slab into the asthenosphere, the positive buoyancy of continental will eventually halt subduction. Already in the early days of plate tectonic theory it has been proposed that when subduction has been brought to standstill, stresses in the subducted plate ultimately lead to break-off of the oceanic slab (McKenzie, 1969). The process of slab break-off shows from seismic tomographic images that indicate gaps at about 200 km depth in the downgoing slab in several subduction zones of the Mediterranean (Spakman et al., 1988;Wortel andSpakman, 1992, 2000).
Numerical models suggest the detachment process takes place by necking of the slab due to tensile stresses, where localized thinning serves as separation mechanism between continental and oceanic lithospheres, with the subsequent sinking of the oceanic fragment (Gerya et al., 2004;Andrews and Billen, 2009; Van Hunen and Allen, 2011).
There are a few continental collision studies that make use of analogue models, which we will discuss hereafter in more detail. We discriminate between models that describe the rheological properties of the applied materials as elasto-plastic, in the sense that these materials should not exhibit viscous behavior, and studies that use visco-elastic materials for modeling lithosphere deformation. Another distinction can be made between studies that use external forcing as driving forces and models that are driven solely by buoyancy forces (Schellart and Strak, 2016).
Analogue models using elasto-plastic materials and external forcing by a piston include studies by Chemenda et al. (1995); Boutelier et al. (2004) and Boutelier and Cruden (2017), where the latter two studies are set up as thermo-mechanical experiments and lead in some of the experiments to break-off of the oceanic slab. As far as we know, the model setup by Boutelier and Cruden (2017) (next to that of Boutelier and Oncken (2011)) is unique in the sense that the horizontal forces exerted by the piston are actually measured, and can thus be checked for proper dynamic scaling. Boutelier and Cruden (2017) describe the applied analogue lithosphere material as elasto-plastic, in other words: there is no description of the material's viscous properties. However the question remains how the strain rate insensitive rheologies applied in these experiments can be related to viscous properties that are commonly attributed to the deeper lithosphere (e.g. Bürgmann and Dresen (2008)). Regard et al. (2003Regard et al. ( , 2005Regard et al. ( , 2008 and Faccenna et al. (2006) use externally forced analogue models to study continental collision using Newtonian viscoelastic materials to represent the lower crust and upper mantle, next to an upper crust that deforms in a brittle manner. In those cases where the slab is allowed to penetrate the lower mantle, the models lead to detachment, in the form of distributed thinning.
The distributed thinning, resulting in dripping of the oceanic lithosphere, may be resulting from both the Newtonian rheology (Schmalholz, 2011), as well as the relatively low viscosity contrast (two orders of magnitude) between the viscous lithosphere and asthenosphere (see for a comparable numerical result Pysklywec et al. (2000)). Guillaume et al. (2013) present a model that includes subduction of continental lithosphere, where the model is driven by buoyancy only and Newtonian viscous materials are applied. However, no break-off occurs after continental material has subducted. In the following we test whether it is possible to create localized slab break-off using the nonlinear viscoelastic material that we have developed.

Model setup
We set up a laboratory model of oceanic and continental subduction, that includes a nonlinear viscoelastic lithosphere and a Newtonian viscous asthenosphere, see figure 10. We present two model setups: one where the model has a long oceanic plate that after subduction becomes supported by the bottom of the model tank (i.e., the upper mantle-lower mantle boundary) before the continent is forced to subduct; in the second setup the model features a very short oceanic slab and a continental slab, designed such that after subduction the slab will hang freely in the surrounding asthenosphere (comparable to the numerical model by Tscharner et al. (2014)), while the bottom of the tank still coincides with the upper to lower mantle transition. We enforce a no-displacement boundary condition at the continental end of the plates for both models. Table  6 provides the dimensions of the models and the length scaling; the chosen length scaling is such that 1 cm in the model represents 33 km in nature. We also ran two comparable models at a smaller length scale (where 1 cm in the model represents 66 km in nature), for which we show the results in the supplement.
We use the organic plasticine-PDMS-silicone oiliron powder mixture 60OPL 40PDMS 24oil to create our oceanic and continental lithospheres, whereas we use glucose with a Newtonian rheology for the low viscous upper mantle. Densities of the oceanic and continental lithospheres are representing bulk lithospheric densities, where we use varying iron powder content to correct densities (see tables ?? and 7). We base the density differences between each lithospheric unit and the asthenosphere on bulk density estimates from Cloos (1993), see equation (11). Table 7 provides all material properties and scale factors between model and nature (using scaling relations from section 3), complemented by realistic natural ranges from relevant literature. Rheological similarity for the model lithosphere (60OPL 40PDMS 24oil) follows from figure 9, which shows that the analogue lithosphere rheologies fall within the estimated range for average lithosphere rheology (section 2.3) for scaling option I, which coincides with the length and density difference scales used here (tables 6 and 7, respectively). Furthermore, table 7 shows that elastic shear moduli µ e for the lithosphere are of the right magnitude compared to natural values; anelastic shear moduli µ a (down to one order of magnitude lower than µ e ) are in the range of estimated natural values (Chopra, 1997). We correct model images for perspective distortion, such that all views are perpendicular to the model axes, and straight lines in the model are indeed depicted as straight lines in the images. From the side view images we extract the side top and bottom edges of the slab using the canny edge detection algorithm (Canny, 1986). From the edge locations we derive the thickness of the slab in time. Model time scales to nature time using equation 12: 1 model minute scales to 0.75 Myr in nature (table  7).

Results of slab break-off model
After we initiate subduction by pulling glucose over the first few cm of the oceanic plate, subduction progresses in both models by rollback of the slab and retreat of the trench. In the long slab model the slab touches the bottom of the tank, after which all negatively buoyant oceanic lithosphere subducts. In this process, a part of the buoyant continental lithosphere is forced to subduct as well (figure 11), and the rollback velocity slows down to zero. Around the passive margin the slab slowly thins, a process that starts to accelerate around 30 minutes (23 Myr scaled) after subduction has started. Figure 13 shows the evolution of the minimum thickness of the slab: thinning starts to accelerate at 28 minutes (21 Myr scaled), followed by initiation of a slab window updip of the passive margin that ultimately leads to full detachment of the oceanic lithosphere around 42 minutes (32 Myr scaled). A movie of the side view images augmented with the detected edges and estimated thickness can be found in supplementary movie S1.
We have performed two models with a short slab (figure 12), that leads to very similar results, with only small differences in the timing of breakoff. Here we discuss only the short slab model 2, as short slab model 1 has less quality images; figures and a movie of the latter can be found in the supplement. In the short slab model the buoyant continental lithosphere is forced to subduct before the oceanic slab touches the bottom of the tank, as the oceanic slab is significantly shorter than the upper mantle depth. Subsequently, around the time that the trench retreat velocity approaches zero (at 20 minutes/15 Myr scaled), the subducted part of the lithosphere is fully suspended from the continental lithosphere at the surface of the model. After 30 minutes (23 Myr scaled), a tear initiates at the right slab edge, just above the oceanic -continental transition, at an equivalent depth of 250 km. While slab break-off commences with the edge tear, the break-off process is complemented by a slab window and wholesale thinning at the same depth. Figure 13 shows accelerated thinning of the slab after 45 minutes (34 Myr scaled), that coincides with the tear migration visible from the front view in figure 12. The total break-off period is 20 minutes (from 30 to 50 minutes model time), which represents 15 Myr in nature.
We can observe that in the long slab model and the two short slab models the final tear occurred in the continental domain, just updip of the oceanic -continental transition. This is also visible in Sup-  (w x l). The thickness of the continental lithosphere tends to reduce after final assembly of the model due to an excess in gravitational potential energy. The thickness of the continental lithosphere may therefore reduce by 25% before the beginning of the experiment. The room temperature during the experiment with the long slab model and short slab model 1 and 2 was 24, 20 and 23 • C, respectively. long slab model short slab model 1 / 2 range in nature ln/lm 160·10 −3 528 · 10 3 160·10 −3 528 · 10 3 l ol 230·10 −3 759 · 10 3 80 / 70 ·10 −3 264 / 231 · 10 3 plement S1 and S2 that display the extracted edge and passive margin location in side view, together with the evolution of the thickness with depth. After break-off, in all models the continental lithosphere rebounds to the surface, while the oceanic lithosphere subducts further until it rests on the bottom of the model upper mantle. The Reynolds number is on the order of 10 −4 , which justifies the assumption that inertia can be neglected. The two long slab models that we have performed at a two times smaller length scale lead to similar subduction and slab rollback, but did not result in any break-off (supplement).

Discussion of the slab break-off model
The location of tearing in both long and short slab models is just above where one would expect the largest tensile stresses for a stationary hanging slab: the transition between buoyant and negatively buoyant lithosphere (e.g. Van Hunen and Allen (2011); Tscharner et al. (2014)). The tearing localizes updip of the continental-oceanic transition while the sinking oceanic plate behaves plate-like.
In other words, the plate does not drip-off such as Regard et al. (2003) observe in analogue models of continental collision using a lithosphere with a Newtonian rheology. Schmalholz (2011) present an analytical study of slab detachment that shows the influence of power-law slab rheology on the minimum slab thickness in time. The acceleration in thinning that we see in figure 13 is very similar to the results of the latter, where such an acceleration in thinning is typical for lithosphere with a power-law rheology, in contrast to linear thinning for a Newtonian slab rheology. We thus account the localization of deformation to the power-law rheology of the lithosphere materials.
For the short slab models it is straightforward to approximate the tensile stresses at the continentaloceanic transition at the time before thinning commences, by calculating the buoyancy forces resulting from the hanging oceanic slab. As the rollback has ended we can assume viscous forces exerted by mantle flow can be neglected. The stress at this location and time can be approximated by: σ = g l ol ∆ρ ol−a which is about 136 Pa. From figure 6 it shows that at this stress the material rheology is in the nonlinear regime, and resulting viscosity is about 1 · 10 6 Pa s. The absence of break-off in models with at a twice smaller scale (supplement) shows the in- Table 7: Model parameters. Density Reference densities of the asthenosphere, and continental and oceanic lithosphere are based on Cloos (1993), where the lithospheric densities represent bulk density. The oceanic lithospheric density of 3310 kg m −3 (∆ρ ol−a = 80 kg m −3 ) assumes an oceanic crust that has been transformed into eclogite during subduction. Values within brackets denote values likely for lithosphere around a continental margin. • To arrive at these densities, 37% and 92% mass fraction iron powder (3.6 and 8.7 volume % respectively) has been used for continental and oceanic lithospheric materials (with 60% organic plasticine, 40% PDMS and 24% silicone oil in terms of mass). Asthenosphere viscosity † Asthenosphere steady-state viscosities estimates for subduction zones are commonly on the order of 10 19 Pa s, based on observations of postseismic deformation (generally estimates are in the range [10 18 -10 20 ] see table 1 of Wang (2007)). This viscosity estimate is in line with results from recent studies that interpret postseismic deformation after recent megathrust earthquakes, such as the 2004 Sumatra-Andaman earthquake (Hoechner et al., 2011;Broerse et al., 2015), the 2010 Maule earthquake (Klein et al., 2016;Bedford et al., 2016) or 2011 Tohoku earthquake (Hu et al., 2016). However, the asthenosphere viscosities reported for subduction zones may be more applicable to the mantle wedge than the suboceanic asthenosphere Hirth and Kohlstedt (2003), for which we adopt the higher value of ≈ 10 20 Pa s (Hu et al., 2004;Broerse et al., 2015). Analogue material rheology ‡ Viscosity at 20 • C according to manufacturer specifications. Lithosphere rheology • The estimate for the mean zero-stress mean viscosity of the ductile lithosphere η0 l is based on diffusion creep parameters from Hirth and Kohlstedt (2003), see explanation in section 2. * Estimates for power-law exponent n for olivine range from 3 to 4.9 (Karato and Wu, 1993;Korenaga and Karato, 2008). The transition stress σt has a large range, but correlates strongly withη0, see figure 1. High values of σt apply to low values ofη0. Elastic parameters mantle lithosphere shear modulus from PREM (Dziewonski and Anderson, 1981). § Anelastic shear modulus values from Chopra (1997)  6.25 · 10 17 ηa [Pa s] 1.6 · 10 2 1.0 · 10 20 10 18 − 10 20 † η0 l 2.0 · 10 8 1.2 · 10 26 2 · 10 22 -4 · 10 28 • nn/nm (eq. 8) 1 n l [ ] 3.6 3.6 3 -4.9 * σn/σm (eq. 11) 1.58 · 10 6 σt l [Pa] 20 3.1 · 10 7 2 · 10 6 − 3 · 10 9 μe 1.8 · 10 4 2.8 · 10 10 6.7 · 10 10 min(µa) 2.8 · 10 3 4.4 · 10 9 2 · 10 9 − 37 · 10 9 § time tn/tm (eq. 12) 3.95 · 10 11 (1 minute model : 0.75 Myr nature)   We estimate the thickness of the slab from the side view images, for the location of top and bottom slab edges and the evolution of the lithospheric thickness we refer to movie S1 and S2 in the supplement. The movie shows that the thinnest part of the slab is generally just updip of the passive margin. As in both models the tear manifests as a slab window that increases in size, the final thickness in these figures represents the width of a thin string of lithospheric material that connects the part of the lithosphere at the bottom of the model and the part that rises back to the surface. In the images the edges appear about four pixels wide, which translates into an uncertainty of located edges of about 0.06 cm. The colored bars denote the epochs of stalling of rollback, initiation of tearing, rapid slab window propagation, and final detachment of the oceanic lithosphere.
fluence of the non-linear rheology on model outcomes: as stress scales with length scale (equation 11) and the lithosphere rheology has a power-law dependence on stress (for the stress levels where break-off occurs), a reduction in scale translates to a power-law reduction in strain rates. The difference in model behavior at varying length scales reiterates that in the design of analogue experiments, the dynamic scaling of a power-law material should be carefully evaluated.
Recently, Boutelier and Cruden (2017) published results on a thermal-mechanical laboratory model of slab break-off, but as they describe the model lithospheric materials as exhibiting a elasto-plastic rheology (i.e., yielding at a specific stress, without viscosity) the results are difficult to interpret knowing that the lithosphere likely has a viscosity. Our results however show that localization can occur in a properly dynamically scaled analogue model using a power-law material, in a tectonic setting without adverse boundary effects such as in previous studies of the influence of rheological nonlinearity on shear localization (e.g. Schrank et al. (2008)).

Conclusions
We have created new materials to represent the lithospheric mantle in analogue models, based on silicone polymers and organic plasticine using creep and recovery tests.
Combinations of organic plasticine and PDMS exhibit deformation rates with a nonlinear dependence on stress. If mixed in a mass ratio 60:40, the plasticine-PDMS mixture behaves as a power-law viscoelastic material with a stress exponent of around 3.5, which fits in the range of reported values for olivine aggregates. At low stresses we observe that the material shows Newtonian behavior, with a transition to power-law viscosity beyond a critical stress level, comparable to composite flow laws for olivine creep that include both Newtonian diffusion and nonlinear dislocation creep. Addition of low-viscosity silicone oil is effective in lowering the material viscosities to suitable values that allow the use of the material in scale models of the ductile lithosphere.
Application of the developed nonlinear material in a model of forced subduction of continental lithosphere, leads to localized deformation and ultimately results in slab break-off.

Acknowledgements
This manuscript benefitted from the accurate comments of two anonymous reviewers. Taco Broerse acknowledges financial support from the ISES program. We greatly appreciate the help of Janet Zulauf who has shared rheological data of organic plasticine and sent samples of plastiline. We would like to thank Benjamin Guillaume as well as João Duarte for sharing information on the use of iron powder in analogue materials. We would like to thank Nicolai Nijholt for comments on a previous version of the manuscript. Antoine Auzemery and Job Arts have been very helpful in the assembly of the laboratory models.

Data availability
The rheometry data have been published in an accompanying data publication: Broerse, Taco

A Strain plots
This section contains the loading and relaxation strain curves for all materials used at the stress levels that have been used in the estimation of rheological parameters.  Figure 15: Strain during loading and relaxation for 60 OPL 40 PDMS at each stress level, with elastic strain γe, anelastic strain γa and permanent strain γp. Permanent strain is used to compute the mean strain rate and steady-state viscosityη (equation (20)).  Figure 16: Strain during loading and relaxation for 60 OPL 40 PDMS 8 oil at each stress level, with elastic strain γe, anelastic strain γa and permanent strain γp. Permanent strain is used to compute the mean strain rate and steady-state viscosityη (equation (20)).  Figure 20: Strain during loading and relaxation for 60 OPL 40 PDMS 24 oil @ T = 30 • C at each stress level, with elastic strain γe, anelastic strain γa and permanent strain γp. Permanent strain is used to compute the mean strain rate and steady-state viscosityη (equation (20)).   Figure 21: Strain during loading and relaxation for 60 OPL 40 PDMS 24 oil 37 Fe at each stress level, with elastic strain γe, anelastic strain γa and permanent strain γp. Permanent strain is used to compute the mean strain rate and steady-state viscosityη (equation (20)).   Figure 22: Strain during loading and relaxation for 60 OPL 40 PDMS 24 oil 92 Fe at each stress level, with elastic strain γe, anelastic strain γa and permanent strain γp. Permanent strain is used to compute the mean strain rate and steady-state viscosityη (equation (20)). In figure 23 we show how we derive the average rheology using emperical parameters for dry olivine and an oceanic temperature depth profile. Here we provide the parameters used, taken from Hirth and Kohlstedt (2003). Figure 23 shows average rheologies for four cases: dry vs. wet olivine rheologies, and continental geothermal gradients.   Hirth and Kohlstedt (2003). We apply a correction to C for parameters based on uni-axial experiments: C = C 2 3 (n+1)/2 (Ranalli, 1995). Grain size d is in µm, stress σ as well as pressure P are in Pa. Temperature T is in degrees Kelvin. Gas constant R is 8.314. We calculate lithostatic pressure P using the values listed for ρ and g. Figure 23: Average lithosphere rheology using wet and dry olivine diffusion and dislocation creep rheological parameters from Hirth and Kohlstedt (2003) applied to continental and oceanic lithospheres (see section 2), left panels: η0 for different values of grain size d, dotted lines denote average η0; right panels: σt, dotted lines denote values below the LAB defined by η0 = 10 20 Pa s. The top depth represents the brittle-plastic transition defined by the intersection of the Goetze criterion and the ductile stress at a strain rate of 10 −14 s −1 .