3-D thermo-mechanical laboratory modeling of plate-tectonics: modeling scheme, technique and ﬁrst experiments

. We present an experimental apparatus for 3-D thermo-mechanical analogue modeling of plate tectonic processes such as oceanic and continental subductions, arc-continent or continental collisions. The model lithosphere, made of temperature-sensitive elasto-plastic analogue materials with strain softening, is submitted to a constant temperature gradient causing a strength reduction with depth in each layer. The surface temperature is imposed using infrared emitters, which allows maintaining an unobstructed view of the model surface and the use of a high resolution optical strain monitoring technique (Particle Imaging Velocimetry). Subduction experiments illustrate how the stress conditions on the interplate zone can be estimated using a force sensor attached to the back of the upper plate and adjusted via the density and strength of the subducting lithosphere or the lubrication of the plate boundary. The ﬁrst experimental results reveal the potential of the experimental set-up to investigate the three-dimensional solid-mechanics interactions of litho-spheric plates in multiple natural situations.


Introduction
Plate tectonic processes are characterized by very large spatial and temporal scales.Consequently, geological data often provide partial insights into their mechanics, and geodynamic modeling, using either experimental or numerical techniques, is routinely employed to better understand their development in space and time.The experimental modeling technique is particularly efficient to investigate threedimensional phenomena (Davy and Cobbold, 1991;Bellahsen et al., 2003;Funiciello et al., 2003;Schellart et al., 2003;Cruden et al., 2006;Luth et al., 2010).However, in multi-Correspondence to: D. Boutelier (david.boutelier@monash.edu)ple experimental models, the rheological stratification of the lithosphere is simplified and the strength variations induced by the temperature gradient through the lithosphere are simulated using various analogue materials with different physical properties (Davy and Cobbold, 1991;Schellart et al., 2003;Cruden et al., 2006;Luth et al., 2010).A drawback of this simplification is that the mechanical properties are retained throughout the entire experiment regardless of temperature variations associated with vertical displacement.
We developed a new apparatus allowing the implementation of three-dimensional thermo-mechanical lithosphericscale models (Fig. 1).In this study, we describe the new temperature-sensitive analogue materials, the advantages of the experimental set-up and the necessary assumptions.We then present experiments illustrating two major features of the apparatus.The first one is the incorporation of a force sensor recording the horizontal convergence-parallel tension or compression in the model lithosphere (Fig. 1) and allowing identification of some key parameters that control the stress regime in the arc and back-arc area during oceanic subduction.The second major feature is that the temperature at D.   3-D sketch of the experimental set-up.Two lithospheric plates made of hydrocarbon compositional systems rest on the asthenosphere modeled by water.A temperature gradient is imposed in the model lithosphere.A thermo-regulator receives the temperature of the water and at the model surface measured by 2 thermal probes and automatically adjusts the length of the heat pulses produced by the lower electric heater and 4 infrared emitters.Plate convergence is imposed at a constant rate and a force sensor installed in the back of the upper plate measures the force in the horizontal convergence parallel direction.Model strain is monitored using a Particle Imaging Velocimetry system imaging the model surface.A second optional camera is employed to follow the model evolution from the side.
the surface of the model lithosphere is controlled and maintained using infrared emitters, which permit an unobstructed view of the model surface (Fig. 1).In turns, this allows using the Particle Imaging Velocimetry (PIV) technique to precisely monitor, from above, the model deformation in time and space (Hampel, 2004;Adam et al., 2005;Hoth et al., 2006Hoth et al., , 2007Hoth et al., , 2008)).

General modeling scheme
The lithosphere is the superficial shell of the Earth capable of undergoing large quasi-rigid horizontal displacements with strain-rates far lower than those experienced by the underlying asthenosphere (Anderson, 1995).This definition provides the general framework for our modeling.Since the viscosity of the asthenosphere is several orders of magnitude lower than the effective viscosity of the lithosphere (Mitrovica and Forte, 2004;James et al., 2009), it can only exert a small shear traction on the base of the lithosphere (Bokelmann and Silver, 2002;Bird et al., 2008).Hence, when studying deformation in a limited area around the subduction zone we can neglect the local effect of viscous coupling between the lithosphere and asthenosphere and replace its larger-scale effect by the action of a piston.Consequently the asthenosphere can be modeled with a low-viscosity fluid (water), whose unique role is to provide hydrostatic equilibrium below the lithosphere.We acknowledge that the role of the asthenosphere is currently simplified (Bonnardot et al., 2008a).The viscous interaction between the lithosphere and asthenosphere will be implemented in a later developmental stage of our experimental set-up.However, before implementing viscous interactions, it is necessary to develop a three-dimensional apparatus standing on the earlier twodimensional thermo-modeling apparatus.
If it is clear that the asthenosphere can be modeled with a fluid, the mechanical behavior of the lithosphere is more complicated.Laboratory measurements of rock strength extrapolated to the conditions of pressure, temperature and strain-rates characteristic of plate tectonics led to the development of the Brace-Goetze strength profile for the rheological stratification of the lithosphere (Goetze and Evans, 1979;Brace and Kohlstedt, 1980;Evans and Kohlstedt, 1995;Kohlstedt et al., 1995).The mechanical behavior of the lithosphere is brittle near the surface and is mainly controlled by frictional sliding (Byerlee, 1978).However at greater depth the lithosphere becomes more ductile, and with further increase of temperature and pressure with depth, the lithospheric behavior becomes more viscous (Ranalli and Murphy, 1987;Ranalli, 1997).One could then represent the oceanic lithosphere with a 3-layer model in which the uppermost layer is elasto-plastic with high brittleness, the second layer is more ductile and the bottom layer is elasto-viscoplastic.In this study, the oceanic lithosphere is simplified and represented by one unique elasto-plastic ductile layer whose strength, however, decreases with depth because of the imposed temperature distribution.The latter is simplified because we impose a linear conductive gradient through the lithosphere (i.e.no heat production), and an isothermal sub-lithospheric mantle (i.e.vigorous convection).
Finally, since we are interested in the large deformation of the lithosphere, we must also consider the material strainsoftening behavior required in order to produce and maintain lithospheric-scale shear zones.In the brittle layers of the lithosphere, localization is a natural consequence of the deformation mechanism.However, strain localization in the ductile regime is not fully understood.It may be promoted by strain softening mechanisms such as dynamic recrystallization and shear heating (Poirier, 1980;White et al., 1980;Rutter and Brodie, 1988;Montési and Zuber, 2002;Hartz and Podladchikov, 2008), but the softening functions associated with these mechanisms are still being debated.An equivalent efficient strain localization process (up to ∼50 % softening) is obtained in our model lithosphere using specially designed strain-softening elasto-plastic materials.

Scaling requirements
The principle of analogue modeling is based on scaling laws, in which geometric, kinematic and dynamic similarity between the natural prototype and the model are achieved by scaling-down lengths, time and forces (Buckingham, 1914;Ramberg, 1967;Davy and Cobbold, 1991;Shemenda, 1994).When scaling relationships are fulfilled, the model behaves mechanically like the prototype but is manageable in a laboratory, and the effects of various key parameters can then be tested using multiple, carefully designed, built and monitored experiments.
For conveniency the model lithosphere must be relatively small.Hence, we decided that 1 cm in the model represents 35 km in nature, which yields a length scaling factor The scaling of the all densities in the model is already set by our choice of using water to represent the asthenosphere.The scaling factor for densities is ρ * = ρ m /ρ n = 3.08 × 10 −1 as 1000 kg m −3 represents 3250 kg m −3 .Consequently, an old, dense, oceanic lithosphere can be modeled using a material with a density of 1040 kg m −3 representing 3380 kg m −3 in nature (i.e. with a negative buoyancy of 130 kg m −3 ), while the low-density continental crust can be modeled with material having a density of 860 kg m −3 representing 2795 kg m −3 in nature.
The scaling of stress is already set by the scaling of hydrostatic pressure ρ g z where depth z scales with length.Since the experiments are produced with normal gravitational acceleration (i.e.g * = 1 or g m = g n = 9.81 m s −2 ), σ * = σ m /σ n = ρ * × L * = 8.79 × 10 −8 and a flow stress of ∼ 10 MPa at the bottom of the lithosphere must be ∼ 1 Pa in the model, while a flow stress of 500 to 1000 MPa in the stronger part of the lithosphere must be 44 to 88 Pa in the model.Therefore, the analogue material employed to model the oceanic lithosphere should have a strength ∼ 1 to ∼ 100 Pa from the bottom to the top.
Before plastic failure, the lithosphere deforms elastically with a shear modulus G n of ∼ 1 to 10 × 10 10 Pa (Dziewonski and Anderson, 1981).Since the dimension of G is that of a stress, it must be scaled by the same ratio: σ * .Therefore, the model shear modulus should be of the order of ∼ 1 to 10 × 10 3 Pa.However, measuring the shear modulus of a very weak material proved to be challenging and therefore the shear modulus could only be measured for low temperatures corresponding to the surface of the model lithosphere.Hence, we must acknowledge that the elastic properties are only approximately scaled.
The scaling of time is chosen in order to scale the temperature variations associated with deformation.The imposed velocity controls the advection of heat in the model, which must be balanced with the diffusion.In order to maintain the same balance between advection and diffusion as in nature, the dimensionless ratio VL/κ, with V being the velocity, L the length, and κ the thermal diffusivity, must be the same in the model and nature (Chemenda et al., 2000).Since the scaling of length has been set already, the thermal diffusivity of the analogue materials controls the scaling of velocity and therefore of time.The later is simply defined in a kinematic sense using the dimensionless ratio Vt/L where t is the time.The scaling of velocity and time is further detailed in Sect.4.5.

Analogue materials
The employed analogue materials are thixotropic dispersions of solid hydrocarbons and powders in oil possessing elastovisco-plastic properties that are strongly dependent on temperature.Like real rocks, the analogue materials exhibit various rheological behaviors for various temperatures and strain-rates, from linear viscous to non-linear viscous to plastic and brittle dilatant properties.The temperature and strain rates of the model lithosphere are chosen so that the materials can be considered purely elasto-plastic (i.e., not sensitive to strain rate) with softening and thus satisfy the scaling requirements presented above.

Composition
The materials are composed of various quantities of paraffin and microcrystalline waxes, Vaseline, and paraffin oil held together with a highly branched alphaolefin polymer.

Mechanical testing
The mechanical properties of the materials have been measured at the Helmholtz-Zentrum Potsdam (GFZ) using a Rheotec RC20 CPS rheometer with cone-and-plate geometry.The materials are liquid at a temperature greater than or equal to the melting point (45-50 • C).To measure the mechanical properties of the materials at temperatures within the range 37-43 • C used in the experiments, the sample is first melted on the heating bottom plate of the rheometer and the rheometer head is slowly lowered into the sample.The temperature is then slowly lowered to 20 • C and the sample is kept at ambient temperature until full strength is developed.The temperature is finally slowly increased to and maintained at the desired temperature of measurement.This protocol allows having full coupling between the sample and both the bottom plate and rheometer head.
The material mechanical behavior at a specific temperature is characterized in simple shear using a series of successive creep tests.A constant shear stress is imposed on the sample for 120 s and then it is ramped up.With low shear stresses, the materials behave elastically and shear strain increases with shear stress (Fig. 2).When the yield stress is reached, shear strain rapidly increases while the shear stress A constant shear stress is imposed on the sample during each time step, and then ramped up.For low shear stresses, the shear strain increases synchronously with the shear stress increase and then remains stable.When the plastic yield stress has been reached, the shear strain increases rapidly during the constant-stress creep step. is maintained constant (Fig. 2).This type of test allows characterizing the elastic shear modulus quantitatively, the plastic yield stress, and the softening/hardening behavior qualitatively .

Elastic properties
The hydrocarbon systems behave elastically for low shear stresses and the shear strain increases instantaneously each time the shear stress is incremented (Fig. 2).Therefore, using the creep tests before failure only, a linear regression of the stress-strain curve provides the elastic shear modulus G, (τ = G × γ ).However, to be measurable, the shear strain increase must be sufficiently large.Large shear stress steps (i.e.≥∼ 5 Pa) are therefore more suitable for quantifying the elastic properties but if the shear stress is increased in too large steps, the precision on the yield stress is not satisfactory.Finally, another restriction on the measurement of the elastic properties arises when the plastic yield stress is small.Then only a few data points can be collected, which is insufficient to derive a meaningful linear regression and thus elastic shear modulus.Consequently, we adopted the following strategy: -For high temperatures (≥ 39 • C), when the material is weak, the mechanical test is performed using small shear stress steps (1 or 2 Pa) and the elastic properties are not derived.
-For low temperatures (≤ 39 • C) the tests are conducted with larger shear stress steps allowing measurement of the elastic properties.
The results reveal that the shear modulus decreases with temperature increase but remains in the range of 1000 to 3000 Pa, which scales to 1.14 to 3.4×10 10 Pa in nature (Fig. 3).Characterisation of the elastic properties thus remains incomplete; however, the presented values of the shear modulus are robust since oscillatory tests (see Boutelier et al., 2008, for precisions on oscillatory tests) also showed a dynamic storage modulus G (equivalent to G in dynamic tests) with the same magnitude (∼ 1500 Pa).

Plastic yield stress
The plastic yield stress is measured using the same series of successive creep tests used to characterize the elastic shear modulus.When the imposed shear stress reaches or overcomes the plastic yield stress, deformation becomes rapidly much larger and it continues during the entire creep step, allowing the rheometer to measure a shear rate.
We qualitatively characterized the softening/hardening behavior using the creep test, during which failure occurred.In this creep test, the shear stress is constant and is equal to the yield stress.We observed that the shear strain increase upon failure is best fitted with a degree-3 polynomial function (i.e.constant acceleration) which is characteristic of plastic softening (Shemenda, 1992).This allows us to affirm that our materials are softening and comparable to the materials developed by A. Chemenda employed in earlier 2-D thermo-mechanical experiments (Chemenda et al., 2000;Boutelier et al., 2002Boutelier et al., , 2003Boutelier et al., , 2004;;Boutelier and Chemenda, 2008).These materials have been shown to soften by 20 to 50 % depending on temperature (Shemenda, 1992(Shemenda, , 1993(Shemenda, , 1994;;Boutelier, 2004).
The rheological stratification of the model lithosphere prior to deformation is obtained when the strength of all the materials constituting the various layers of the lithosphere is known for all the temperatures imposed inside the model lithosphere.In the presented experiments, the lithosphere is made of one single mantle layer and therefore the strength of one single material must be known.Figure 4 shows the variations of the plastic yield stress of our material for the temperatures 38 to 43 • C.However, in the presented experiments the surface temperature was 39 • C and the temperature of the asthenosphere was 42 • C. Within this temperature range, the strength decreases from 50 to 6 Pa (see inset in Fig. 4), which corresponds to a decrease from 5.7 × 10 8 to 6.8 × 10 7 Pa in nature.
Our materials exhibit ductile plastic failure, as permanent deformation can be observed prior to failure in the stresstime diagram and ductile flow can be observed around the fault zone.However, the amount of plastic deformation prior to failure decreases with the temperature.Therefore, our materials become more brittle (i.e.deformation is more localized) with lower temperature.Brittle plastic failure should prevail near the model surface.However, it is not presently possible to implement brittle behavior in our experiments because brittleness is only obtained for very low temperatures, when the materials are too strong.Furthermore, even if we were able to modify our materials in order obtain more brittle behavior near the surface together with a low strength, we would not be able to reproduce the strength increase with depth in this part of the model as our materials are sensitive to temperature but not pressure.We acknowledge that the strength envelope is approximated near the model surface, and that deformation of the near surface is best modeled with granular materials (e.g.Corti, 2008;Malavieille and Trullenque, 2009;Malavieille, 2010).

Thermal diffusivity
The determination of the material's thermal diffusivity is fundamental since it conditioned the scaling of rate and time in the experiments.This parameter is measured using a 1-D cooling approximation.A large sample is brought to a high temperature above the melting point and allowed to cool down.The heat loss is restricted to the upper surface only and the temperature is monitored at a known depth below the centre of the upper surface.In these conditions, a 1-D halfspace cooling approximation is reasonable and the data can be fitted with the analytical solution (Turcotte and Schubert, 1982): where T (z,t) is the temperature at depth z and time t, T e is the external temperature, T i is the initial temperature and κ is the thermal diffusivity.For our compounds, the best fit is obtained with a thermal diffusivity of 2.8 × 10 −8 m 2 s −1 (Fig. 5).Assuming that the modeled rocks have a thermal diffusivity of 1×10 −6 m 2 s −1 (Turcotte and Schubert, 1982), then κ * = 2.8 × 10 −2 and a scaling factor for velocity .8 × 10 4 can be derived from the dimensionless ratio VL/κ.This scaling factor means that a natural subduction velocity of 8 cm yr −1 (i.e.2.54 × 10 −9 m s −1 ) is 0.25 mm s −1 in the model.The scaling factor for time is therefore t * = t m /t n = L * /V * = 2.92×10 −12 , which implies that 1 Ma in nature (i.e.3.15 × 10 13 s) is 92 s in the model.

Experimental apparatus
The experimental set-up comprises a polycarbonate tank 50 × 50 × 30 cm 3 filled with water representing the asthenosphere (Fig. 1).The model lithospheric plates are built in a mold 40 × 40 × 6 cm 3 .The experimental tank is sufficiently larger than the model plates (5 cm on each side) that the sides can be considered to be free.The experimental set-up includes two new features that are presented here.The surface temperature is maintained without obstructing the view of the model surface for optical strain monitoring and the horizontal convergence parallel force is measured at the back of the upper plate.

Particle imaging velocimetry and surface heating
A major difficulty with 3-D experimental models is monitoring the model deformation in its center.Because we wish to explore three-dimensional processes, the model is likely to include along-strike variations of the initial conditions which make the deformation process different in the center and along the sides.Therefore the model side views are not sufficient and it is necessary to monitor strain from the top.Precise spatio-temporal strain monitoring is obtained using the Particle imaging velocimetry (PIV) technique (Hampel, 2004;Adam et al., 2005), a non-intrusive method for accurate measurement of instantaneous velocity/displacement fields using an image correlation technique.Our PIV system is equipped with 10 megapixel cameras enabling a spatial resolution of the displacement below 0.1 mm while the temporal resolution is 0.1 s.However, in the presented experiments the successive PIV images are taken at a time interval of 2 to 5 s, which is sufficient to monitor the slow model deformation (< 0.25 mm s −1 ).
To obtain satisfactory measurements of the displacement, the PIV system requires an unobstructed view of the deforming surface and this surface must have a specific pattern allowing small image sub-samples to be shifted incrementally in the x-and y-directions and correlated between two successive images.It is the correlation of the specific pattern that yields the local displacement vector averaged over the surface-area of the image sub-sample (Adam et al., 2005).To obtain this specific pattern suitable for the correlation technique, dark particles have been sifted on the model surface.The difficulty with the use of the PIV technique with thermo-mechanical models was to impose the surface temperature T s while maintaining an unobstructed view of the model surface.We solved this problem using infrared emitters coupled to a thermal probe and a thermo-regulator (Fig. 1).Four 250 V/250 W infrared emitters equipped with large diffusers were placed 60 cm above the 4 corners of the experimental tank and are oriented towards the center of the model surface.The infrared emitters do not produce any visible light and the PIV cameras do not see the emitted infrared waves.Therefore, the heating system is perfectly invisible to the PIV strain monitoring system.The surface temperature is continuously measured by a thermal probe in one location and the temperature value is given to the thermo-regulator which, depending on the temperature difference between the set temperature and measured temperature, adjusts the length of the pulses emitted by the infrared bulbs.This set-up allows having a surface temperature field that is relatively constant (±0.1 • C) in time and homogeneous (±0.2 • C) in space (Fig. 6).A similar system controls the temperature of the model asthenosphere, however the heating element is a simple 250 V/2000 W electric resistance placed at the bottom of the tank, as in previous thermo-mechanical experimental set-ups (Chemenda et al., 2000;Boutelier et al., 2002Boutelier et al., , 2003Boutelier et al., , 2004;;Boutelier and Chemenda, 2008).

Force monitoring
In order to better understand the stress regime in the arc area during the processes of subduction and/or arc-continent collision, we placed a force sensor in the back of the up-per plate.The upper plate rests against a vertical plate attached to the back-wall of the experimental tank via a 2.5 N force sensor.This set-up allows measuring the force in the horizontal convergence-parallel direction.In this study we present a quantitative analysis of the effects of slab buoyancy, flexural rigidity and interplate friction on the interplate stresses and produced stress/strain regimes in the arc area during oceanic subduction.However, the measurement of the horizontal convergence-parallel stress is only really useful in a quantitative manner when the model is cylindrical and the same process occurs at the same time across the entire width of the model.It is therefore best used in preliminary experiments such as presented in this study.The rationale for placing this force measurement system is that the effects of some parameter values (densities, temperatures) on the interplate stresses can be measured in cylindrical experiments and then the same effects can be assumed in more complex non-cylindrical models.
Our modeling apparatus is computer-controlled and a single computer simultaneously reads the force measured at the back of the upper plate and imposes the rate of the piston.The force sensor is interrogated 30 times per second and the computer can either store the values and maintain convergence rate at a set value or it can adjust the convergence rate in order to achieve a set force.We have not yet fully tested this system and there are parameters that remain to be investigated/tuned.For example, we presently don't know the maximum acceleration that should be allowed or how we should filter the force signal in real time in order to ignore the highfrequency force variations that do not pertain to the model mechanics.We will further discuss the future development of the modeling apparatus and in particular of dynamic models in the discussion section.

Simple oceanic subduction experiments
Intra-oceanic subduction experiments comprise two model lithospheric plates, each made of one single mantle layer.The overriding plate is 20 × 40 × 2 cm while the subducting plate is 25 × 40 × 2 cm.The plate boundary position and geometry are imposed and convergence is orthogonal to the direction of the trench.The experiments are cylindrical and can thus be considered quasi-two-dimensional despite the relatively large model width.Simple cylindrical intra-oceanic experiments are performed for two principal reasons.First, the experiments are produced because the resulting scenarios are more predictable and therefore these experiments allow testing the stress and strain monitoring capabilities of the experimental set-up.The second reason is that the normal stress measured at the back of the overriding plate can be more directly related to the stress conditions (σ and τ ) exerted on the plate boundary, which depends on several parameters such www.solid-earth.net/2/35/2011/Solid Earth, 2, 35-51, 2011 as the density and flexural rigidity of the lower plate or the interplate friction.Here we present two simple experiments in which we varied the density of the subducting slab and the interplate friction.
In Experiment 1, the subducting lithophere is neutrally buoyant (ρ l = ρ a ) and the plate boundary is lubricated with a 0.1 mm-thick layer of parafine oil (η = 0.5 Pa s).We assume that in this situation the shear traction is null.Figure 7 shows the model surface photographs at various stages with the velocity vectors derived from the image correlation technique and the incremental (i.e. between 2 successive images) convergence-parallel normal strain (E xx ).The PIV monitoring yields vectors that are very honogeneous in both direction and magnitude within the subducting plate (Fig. 7a-c).The model plate therefore moves as a quasi-rigid body.The upper plate does not move and is only slightly deformed during subduction initiation (Fig. 7a).Therefore plate convergence is only accomodated by sliding along the interplate zone.In the conditions of this experiment, the normal stress measured at the back of the upper plate is the horizontal compression due to the flexural rigidity of the plate (Shemenda, 1993).During subduction initiation, the stress σ xx increases to a value of ∼ 12 Pa and then remains approximately at this level during the rest of the experiment (Fig. 8).This normal stress is due to the non-hydrostatic normal stress σ n exerted on the plate boundary by the subducting plate because it resists bending.Knowing the value of the normal stress σ xx within the upper plate we can estimate the horizontal component (F p ) h of the pressure force F p .Since we also know the inclination α, width W and thickness H of the plate, we can estimate the magnitude of the pressure force F p and the depth-averaged non-hydrostatic normal stress σn from which it derives (Fig. 10a): Note that it is possible to increase or decrease the flexural rigidity of the subducting plate and thus the depth-averaged non-hydrostatic normal stress by adjusting, for example, the temperatures (and thus strength) in the model.
In Experiment 2, the subducting lithosphere is denser than the asthensphere by 4 %.Consequently, the slab exerts a pull on the plate boundary which subtracts from the compression due the flexural rigidity.Also, since the viscosity of our model asthenosphere is very low and the trailing edge of the plate is fixed, the negative buoyancy of the subducting lithosphere leads to the formation of a very steep slab (Fig. 9).In this experiment, we also increased the interplate friction by sifting a small amount of sand grains at the surface of the subducting plate.Therefore, the normal stress σ xx recorded at the back of the upper plate reflects the sum of the compression due the flexural rigidity, the tension due to the slab pull and the compression due to the increased friction.Of these 3 key parameters, only one is known from Experiment 1.However, the effects of the slab negative buoyancy and high interplate friction evolve differently in time.At the begining of the experiment, the effect of the slab pull force is minor because the slab is not developped yet.Consequently we recorded a normal stress σ xx that is larger than that recorded in Experiment 1 (Fig. 8).Assuming that the effect of the slab pull force can be neglected at the begining of the experiment, the stress difference between experiments 1 and 2 must be due to increased friction.We can therefore estimate the horizontal component of the friction force F f (Fig. 10c): then and finally the depth-averaged shear stress After a short plateau, the recorded normal stress σ xx decreases and it keeps decreasing until the end of the experiment.The steep geometry of the subducted slab (Fig. 9) reveals that the decrease of the horizontal compression is due to the slab pull force exerted by the subducting slab on the interplate zone.In the presented experiment the horizontal stress becomes very close to zero near the end of the experiment.In other similar experiments it has been possible to obtain a tension.However, the magnitude of the tension cannot be very large because the plates are not strongly attached to the piston and back wall and would, under the effect of a large tension, detach and move towards the center of the experimental tank.σ xx is due to the combined effects of the horizontal component of the pressure force due to the slab negative buoyancy (F p2 ), the horizontal components of the pressure force due to flexural rigidity (F p1 ) and the friction force (F f ): We estimate the pressure force due to the slab negative buoyancy F p2 (Fig. 10b): and the associated depth-averaged non-hydrostatic normal stress σn σn = Note that for clarity only 1/4 × 1/8 of all the velocity vectors are presented here.From the velocity field the normal strain in the convergence direction E xx is derived.The magnitude of E xx is calculated using the displacement field computed between 2 successive images taken 5 s apart (i.e.divide by 5 to obtain time-averaged strain rate).Both the vector field and the strain field reveal that the both plates are largely undeformed and convergence is accommodated by sliding of the lower plate under the upper plate.The presented data clearly show that (1) the experimental apparatus allows modeling plates which move with little/no internal deformation, (2) the monitoring technique allows precise spatial resolution of model deformation, and (3) the measured stress allows defining 3 end-member subduction regimes characterised respectively by compression due to the bending strength of the subducting plate, tension due to the negative buoyancy of this plate and compression due to the interplate friction (Fig. 10).

Forced subduction initiation experiment
The presented subduction initiation experiment contains only one model lithospheric plate with dimensions 40×40×3 cm.A notch is created in the middle of the plate's length which runs across the entire plate's width.In the notch the plate thickness is only 2 cm.It constitutes a weak zone in the model lithosphere because the thickness is reduced, and the  10.Sketches of the main subduction regimes and the key parameters causing them.A compressive subduction regime is caused by the flexural rigidity of the lower plate which resists bending (a).However, this regime is only obtained when the slab is neutrally buoyant because the negative buoyancy of the subducted lithosphere exerts a tensile non-hydrostatic normal stress on the plate boundary which generates a horizontal tension in the upper plate (b).Finally, a horizontal compression can also be generated because of a large interplate friction (c).thermal gradient is higher, which makes the model lithosphere at shallow depth above the notch weaker than elsewhere in the model.The model plate is shortened at the same constant rate as in Experiments 1 and 2, which results in deformation of the plate near the notch and finally the formation of a subduction zone (Fig. 11).We do not pretend that the presented mechanism is at play in nature during subduction initiation.This experiment is realized in order to observe how strain localization due to plastic strain-softening can lead to the formation of a lithospheric-scale shear zone in our model lithosphere.Furthermore, we use this experiment to illustrate how the PIV system allows having a precise monitoring of when the various deformation structures are active.
The experiment being cylindrical is also monitored from the side.Successive side views revealed that shortening is accommodated by two conjugate shear zones nucleated from the roof of the notch and resulting in the formation of a popup above it (Fig. 11a).With further shortening, the shear zone located left of the notch becomes dominant, accommodates most of the shortening and thus forms a new subduction zone (Fig. 11b-f).However, the side views also reveal that both plates underwent some shortening near both the piston and the back wall (Fig. 11b).
The PIV monitoring of the model surface allows us to describe more precisely the evolution of model deformation in space and time.The shortening indeed started with the formation of two oppositely dipping shear zones around the notch (Fig. 12a).However, before the shear zone on the left side of the notch became dominant, another shear zone is also created near the back wall (Fig. 12a).At this point, none of the shear zones appear to run entirely across the width of the model plate.This shortening is accommodated by the network of shear zones and distributed across the entire model.PIV monitoring reveals that this process of strain accommodation does not perpetuate.The shear zone located to the right of the notch becomes inactive while the shear zones located left of the notch and near the back wall propagate laterally (Fig. 12b).It is the advantage of using the PIV technique to be able tell when, during the experiment, a structure becomes inactive and when another one becomes active or more important.The shortening is accommodated by fewer and fewer structures as strain weakens the material.Finally the shear zone located near the back wall also becomes inactive (Fig. 12c) and shortening is solely accommodated by one shear zone which evolves into a new subduction zone (Fig. 12c-d).The strain localization process is accompanied by a diminution of the stress recorded at the back of the upper plate once the main fault zone runs across the entire model width and has recorded sufficient strain to be significantly weakened (Fig. 13).The recorded diminution is also due to the fact that, since the mantle lithosphere material is slightly denser than the asthenosphere, the formation of a slab leads to a slab pull force whose action diminishes the horizontal compression.

A framework for modeling lithospheric plate interactions in a convergent setting
The modeling framework presented in this study, with its set of simplifications and assumptions, is designed to investigate the three-dimensional thermo-mechanical interactions between lithospheric plates.Targeted geodynamic problems include the lateral (i.e.along trench) propagation of an arc-continent collision (e.g.Taiwan, Timor, Papua New Guinea) or the deformation of the fore-arc/arc/back-arc along a curved plate boundary (e.g.Central Andes).In collision zones, the deformation of the lithosphere is threedimensional and controlled by the stress conditions exerted on the plate boundary (i.e.non-hydrostatic normal stress and shear traction).In order to better understand such complex three-dimensional situations, scaled three-dimensional modeling of the lithospheric plates interaction is required.Furthermore, since the targeted geodynamic problems are likely to include large vertical movements of the crust and/or lithosphere and hence significant temperature variations, it is important that the modeling framework includes the spatial and temporal strength variations associated with temperature changes.
The presented experiments are large but cylindrical models since we did not implement any lateral (i.e.along-strike) variations of the model structure, mechanical properties or any other initial or boundary conditions.We presented cylindrical models because such simpler models allow us to introduce the key parameters controlling the solid-mechanics interaction of the lithospheric plates during the process of oceanic subduction.Furthermore, these experiments permit detailing the advantages and limitations of our approximations and modeling set-up.However, the new apparatus clearly has the potential for modeling the above-mentioned non-cylindrical geodynamic problems.The present study should therefore be seen as a first step in a series designed to study such three-dimensional non-cylindrical processes, rather than the final step.

Quantitative stress and strain monitoring
One key advantage of the employed experimental setup is the ability to precisely monitor the model deformation from above using the 2-D PIV technique.Because we aim at producing complex 3-D models, we cannot rely on side observations of the model as in previous 2-D thermo-mechanical experiments (Boutelier et al., 2002(Boutelier et al., , 2003(Boutelier et al., , 2004)).The high spatial and temporal resolution of the PIV technique allows a more precise monitoring of the model deformation than conventional strain marker analysis (Boutelier and Cruden, 2008).Theoretically, it is also possible to monitor the vertical displacement of the model surface using 2 cameras providing a stereoscopic view of the model surface (Adam et al., 2005;Riller et al., 2010).This new advance in monitoring technology would be very advantageous because the vertical motion of the model surface can be linked to the distribution of stresses along the plate boundary (Shemenda, 1992), compared with numerical modeling results (Hassani et al., 1997;Bonnardot et al., 2008a,b) or compared with natural data such as long-term uplift/subsidence derived from sedimentary record (Matsu'ura et al., 2008(Matsu'ura et al., , 2009;;Stefer et al., 2009;Hartley and Evenstar, 2010), short-term uplift/subsidence derived from GPS/Coral reefs dating (Taylor et al., 2005;Matsu'ura et al., 2008Matsu'ura et al., , 2009) ) or gravity anomalies (Shemenda, 1992;Song and Simons, 2003).However, the produced topography is very small (∼ 1-2 mm) and the distribution of passive markers is not sufficiently large to permit reliable stereoscopic imaging and therefore vertical motions.A possible workaround would be to employ a three-dimensional video laser to scan the model surface at regular time intervals during the experiment and convert the data to digital elevation models (Willingshofer et al., 2005;Willingshofer and Sokoutis, 2009;Cruden et al., 2006;Pysklywec and Cruden, 2004;Luth et al., 2010).The strain monitoring system in its present development allows precise spatio-temporal resolution of the 2-D deformation of the model surface.This allows tracking when, during the experiment, the various deformation structures are active.Coupled with stress monitoring and the ability to produce vertical sections through the frozen model in the end of the experiment (Boutelier et al., 2002(Boutelier et al., , 2004;;Boutelier and Chemenda, 2008), the strain monitoring system provides the required data to understand how stress and strain develop in the model.

Imposing a subduction regime in non-cylindrical thermo-mechanical experiments
The stress monitoring system allows investigating parameters that control the subduction regime in cylindrical experiments and then the same effects can be assumed in more complex, non-cylindrical, three-dimensional models.For example, kinematic models of the Andes show that trenchparallel normal shortening as well as trench-perpendicular  normal shortening existed in the center of the plate boundary curvature while the Andes were built (Kley, 1999;Hindle et al., 2002Hindle et al., , 2005;;Oncken et al., 2006;Arriagada et al., 2008;Gotberg et al., 2009).Two-dimensional, numerical simulations reveal that trench-parallel compression is produced near the symmetry axis of a seaward-concave plate boundary if interplate friction is high and/or the subducting lithosphere has a low flexural rigidity (Boutelier and Oncken, 2010).In contrast, trench-parallel compression is reduced along the oblique parts of the plate boundary.However, both the stress conditions on the interplate zone and the 3-D geometry of this zone control whether the trench-parallel stress in the centre of the curvature is a tension or compression.Low dip angle and high convergence obliquity angle favour trenchparallel compression.In order to investigate how deformation develops in the centre of the curvature, it is necessary to produce 3-D modeling with the proper stress conditions along the plate boundary.In the experiments presented in this study, we have explored some parameters that control the stress conditions along the plate boundary.In experiment 1, the non-hydrostatic normal stress is high (12 Pa equivalent to ∼ 120 MPa in nature) because the subducted slab is neutrally buoyant.To favor trench-parallel normal shortening in the centre of the curvature, we should therefore rather use a denser material to model the subducted lithosphere.Indeed, Experiment 2 showed that the effect of the slab pull force due to the negative buoyancy of the subducted lithosphere was to reduce the non-hydrostatic normal stress.However, the density should not be too high or trench-perpendicular tension would be produced.In this case, the model interplate stress conditions would not correspond to the interplate stress conditions that we hypothesized for the Andes and whose effects we wish to investigate.Finally, we can increase the interplate friction and thus the shear traction acting on the interplate zone.In Experiment 2 the depth-averaged shear traction was 13 Pa (equivalent to ∼ 130 MPa in nature), which is very high and would therefore likely produce deformation of the thinner, hotter and hence weaker arc/back-arc (Currie and Hyndman, 2006;Currie et al., 2008).However, shear traction would have to be imposed on the upper shallow part of the plate boundary where the contact between the plates is seismogenic and thus frictional.It is therefore possible, using simple cylindrical experiments as presented in this study, to tune the value of some parameters such as the density or thickness of the subducted plate, or the interplate friction, in order to produce a specific subduction regime characterized by its stress conditions on the interplate zone.The effect of a specific subduction regime can then be tested in noncylindrical experiments where the 3-D geometry of the plate boundary (i.e.dip angle and convergence obliquity angle) can be changed.

Towards dynamical thermo-mechanical experiments
Our experimental set-up will be further developed and will progressively evolve toward a fully three-dimensional, thermo-mechanical and dynamical experimental set-up.The development path includes several milestones.The first one, presented in this study, was to develop a threedimensional thermo-mechanical set-up based on the previous two-dimensional setups (Chemenda et al., 2000;Boutelier et al., 2002Boutelier et al., , 2003Boutelier et al., , 2004)).The next developmental stage will include proper viscous interaction between the lithosphere and asthenosphere in a kinematic framework and in the last stage we will produce dynamic experiments with a constantforce boundary condition.
Since we want to continue developing the presented modeling framework, we must keep the same density of the asthenosphere (ρ a = 1000 kg m −3 ) and adjust the viscosity of the fluid to a specific value derived from the existing scaling of viscosity and the assumed value of the asthenosphere viscosity in nature.The scaling factor for viscosity η * can be deduced from our scaling factors for stress σ * and for time t * .With the values presented in this study we found that η * = 2 × 10 −19 .Therefore if we assume that the asthenosphere has a viscosity of 1 × 10 19 Pa s, we must use a fluid with a viscosity of 2.57 Pa s, and if the asthenosphere viscosity is assumed to be 1 × 10 21 Pa s, we must use a fluid with a viscosity of 257 Pa s.Using a characteristic length of 2 cm corresponding to the thickness of the oceanic lithosphere in our model, we can can compute the Reynolds numbers (R e = ρV L/η) associated with each end-member value of the asthenosphere viscosity.For a model asthenosphere of 2.57 Pa s corresponding to 1 × 10 19 Pa s in nature, R e = 1.95 × 10 −3 , while for a model asthenosphere of 257 Pa s corresponding to 1 × 10 21 Pa s in nature, R e = 1.95 × 10 −5 .Therefore R e 1, the viscous forces are large compared to inertia and the flow patterns in the mantle due to the motion of the plates should be properly reproduced.
The second developmental step will be to impose a constant-force boundary condition in order to produce dynamical experiments.Our modeling apparatus is computercontrolled and the computer can already either store the force value and maintain a set convergence rate or it can adjust the convergence rate and store the displacement, in order to achieve a set force.This feature will enable us to perform dynamic experiments, taking into account the effects of velocity variations on the thermal structure and thus strengths in the model.
Experiments will be conducted using various boundary conditions.For example, the constant force measured in a oceanic subduction experiment with a neutrally buoyant slab can be imposed in a second run as a constant-force boundary condition.We expect to produce the same subduction velocity during the oceanic subduction stage that produced the proper thermal structure in the subducted slab.Furthermore, if we also introduce a continental passive margin, we expect a velocity reduction during the burial of the buoyant continental crust, whose effects in a three-dimensional thermomechanical set-up are poorly constrained.More complex types of boundary condition used in some recent numerical simulations (i.e.constant velocity for a certain time or amount of convergence, followed by a constant force/stress) (Faccenda et al., 2008(Faccenda et al., , 2009) ) should also be possible.
Implementation of viscous interaction between the lithosphere and asthenosphere is just a matter of choosing a new material that has the appropriate physical properties to achieve complete compliance of the scaling rules and of the force balance.This will only be necessary for future experiments studying a more complete system, but does not pertain to the experimental apparatus and modeling strategy, the topic of this paper.The implementation of the constant force boundary condition will require further development of the experimental setup, but the basic driving system has already been produced.What remains to be found is only a technical solution in order to precisely control the force or rate and minimize fluctuations.This is not a fundamental modification to the setup but only a minor technical improvement.

Conclusions
We have demonstrated, in this study, the potential of a newly developed three-dimensional, thermo-mechanical analogue modeling setup to investigate complex, three-dimensional, geodynamical problems.The targeted situations include the deformation of the fore-arc, arc and back-arc along a seaward-concave plate boundary such as in the central Andes.The deformation near the symmetry axis is fundamentally three-dimensional with both trench-parallel and trenchperpendicular normal shortening (Kley, 1999;Hindle et al., 2002Hindle et al., , 2005;;Arriagada et al., 2008), and is most likely due to the stress conditions along the plate boundary.We have shown that we can estimate these stress conditions in our analogue models using simple cylindrical experiments.Furthermore, varying some parameter values such as bending strength and relative buoyancy of the lower plate or the interplate friction, we can control these stress conditions along the plate boundary and impose them in non-cylindrical models.
We also demonstrated that the deformation resulting from these imposed stress conditions on the interplate zone can be precisely monitored using the PIV system.This strain monitoring system allows characterization and quantification of horizontal deformation (i.e.E xx , E yy , E xy and E yx ), while model sections after deformation provide access to the final vertical deformation (i.e.amount of thinning or thickening in vertical section E zz ).Furthermore, the high spatial and temporal resolution of the PIV system allows tracking the propagation of the deformation.Such feature is particularly useful for investigating arc-continent or continent-continent collisions, which generally initiate in one location and propagate laterally along the plate boundary, as in Taiwan (Suppe, 1984), Timor (Searle and Stevens, 1984;Harris, 2011), or the Urals (Puchkov, 2009).
The modeling framework presented in details in this study and including new temperature-sensitive elasto-plastic analogue material as well as a new modeling apparatus with force monitoring and precise strain monitoring posesses the potential to be the foundation for multiple investigations into the complex 3-D interactions between lithospheric plates.
Fig.1.3-D sketch of the experimental set-up.Two lithospheric plates made of hydrocarbon compositional systems rest on the asthenosphere modeled by water.A temperature gradient is imposed in the model lithosphere.A thermo-regulator receives the temperature of the water and at the model surface measured by 2 thermal probes and automatically adjusts the length of the heat pulses produced by the lower electric heater and 4 infrared emitters.Plate convergence is imposed at a constant rate and a force sensor installed in the back of the upper plate measures the force in the horizontal convergence parallel direction.Model strain is monitored using a Particle Imaging Velocimetry system imaging the model surface.A second optional camera is employed to follow the model evolution from the side.

Fig. 2 .
Fig. 2.Procedure employed to measure the material's mechanical properties.A constant shear stress is imposed on the sample during each time step, and then ramped up.For low shear stresses, the shear strain increases synchronously with the shear stress increase and then remains stable.When the plastic yield stress has been reached, the shear strain increases rapidly during the constant-stress creep step.

Fig. 3 .
Fig.3.Strain-stress curve obtained using the shear strain augmentation synchronous with the shear stress growth prior to reaching the plastic yield stress.

Fig. 4 .
Fig. 4. Evolution of the plastic yield stress with temperature.The plastic yield stress decreases with increasing temperature.Inset:Once the temperature gradient in the model lithosphere is known, the experimental curve is used to draw the strength envelop of the model lithosphere prior to deformation.

Fig. 5 .
Fig. 5.The material thermal diffusivity is fitted with a 1-D cooling solution.The best fit provides a thermal diffusivity of 2.8 × 10 −8 m 2 s −1 .

Fig. 6 .
Fig. 6.Infrared image of the model surface.The homogeneity of the surface temperature T s is measured using an infrared camera.After 3 h of heating, the temperature field shows only one small ridge associated with the plate boundary.The temperature variations across the model surface are otherwise smaller than 0.2 • C. In this test the surface temperature was 38 • C and the temperature of the asthenosphere was 42 • C. In the bottom right darker corner, the model is masked.

Fig. 7 .
Fig. 7.Surface views of the model in Experiment 1. Correlation of successive images allows the derivation of a velocity field (arrows).Note that for clarity only 1/4 × 1/8 of all the velocity vectors are presented here.From the velocity field the normal strain in the convergence direction E xx is derived.The magnitude of E xx is calculated using the displacement field computed between 2 successive images taken 5 s apart (i.e.divide by 5 to obtain time-averaged strain rate).Both the vector field and the strain field reveal that the both plates are largely undeformed and convergence is accommodated by sliding of the lower plate under the upper plate.

Fig. 8 .
Fig.8.Convergence parallel normal stress recorded in the back of the overriding plate in Experiments 1 and 2. In Experiment 1, the recorded compressive stress is due to the flexural rigidity of the plate which resists bending.In Experiment 2, the stress at the beginning of the experiment is due to the combined effects the flexural rigidity and high interplate friction.However, during Experiment 2, the compressive stress decreases because the subducted slab becomes longer and its pull becomes stronger.

Fig. 9 .
Fig. 9. Side images of the model in Experiment 2. The negative buoyancy of the subducted lithosphere tends to steepen the slab.

Fig. 11 .
Fig. 11.Side views of Experiment 3. The model is composed of one single layer of mantle lithosphere with a notch in the middle of the plate's length.The side pushed by the piston (i.e.right-hand side) becomes the upper plate while the other side slides down a newly created subduction zone.

Fig. 12 .
Fig.12.Views of the model surface in Experiment 3, with the velocity vectors and convergence parallel normal strain E xx .The magnitude of strain is calculated for a time increment of 2 s between correlated images.The active faults are drawn with solid line while inactive faults are drawn with a dashed line.The PIV images confirm that the shortening is initially accommodated by a pop-up located above the notch with two opposite verging thrusts rooting near the roof of the notch (a).The PIV images also reveal that the fault located near the back wall is important at the beginning of the experiment (a, b) but becomes inactive while the future main fault develops (c, d).

Fig. 13 .
Fig.13.Convergence parallel normal stress recorded in the back of the overriding plate in Experiment 3. The stress keeps increasing at the beginning of the experiment despite the formation of multiple shear zones.The stress decreases only once one shear zone becomes through-going.Strain is then localized into this zone (see Fig.11b) and stress decreases.

Table 1 .
Parameter values adopted for the models ( a refers to Exp. 2 and b refers to Exp. 3), scaled to nature and scaling factors.The plastic strength of the materials decreases with depth in each layer.In this table we indicate the strengths averaged over the layer thickness.