2-D finite displacements and strain from particle imaging velocimetry (PIV) analysis of tectonic analogue models with TecPIV

2-D finite displacements and strain from particle imaging velocimetry (PIV) analysis of tectonic analogue models with TecPIV

Boutelier, Schrank and Regenauer-Lieb

Solid Earth, 10, 1123–1139, 2019

We published a paper detailing the progress made on TecPIV:  2-D finite displacements and strain from particle imaging velocimetry (PIV) analysis of tectonic analogue models with TecPIV. The open-access paper can be downloaded from the link above.

More information on the software here, and the source here.

TecPIV outputs with strain rate components instead of velocity gradients



Finite deformation ellipses from cumulative Lagrangian sums.

New paper on subduction initiation

Can subduction be initiated at a transform fault? The short answer is probably not. The not so short answers are below and in the linked paper.


Initiation of Subduction Along Oceanic Transform Faults: Insights From Three-Dimensional Analog Modeling Experiments

David Boutelier and David Beckett
School of Environmental and Life sciences, University of Newcastle, Newcastle, NSW, Australia

Subduction initiation is a fundamental component of the plate tectonic theory, yet how subduction starts remains controversial. Oceanic transform faults and fracture zones have been proposed as sites of subduction nucleation because they are thought to be mechanically weak and the large buoyancy gradient across these faults because of the difference in the age of the lithosphere, was thought to facilitate foundering. Self-sustaining subduction, defined as subduction driven by the negative buoyancy of the sinking lithosphere might be achieved if at least ~100 to 150 km of convergence can be imposed on an oceanic fracture zone with sufficient buoyancy gradient across the fault.

Previous modelling however did not take into account the fact that the age of the lithosphere and therefore its strength and buoyancy not only varies across the oceanic fault zone, but also along the strike of the fault since many oceanic transform fault link segments of spreading ridges. Here we investigate using three-dimensional analog models how the spatial distribution of strength and buoyancy along and across an oceanic transform fault zone affects the polarity of subduction and whether self-sustaining subduction can be obtained. We designed three-dimensional analog experiments in which two oceanic lithospheric plates are separated by a weak transform fault and convergence is imposed in the horizontal direction perpendicular to the strike of the fault. The spatial distribution of plate thickness and buoyancy are varied along and across the strike of the transform fault, and whether self-sustaining subduction is obtained is assessed using a force sensor.

Cylindrical experiments reveal that subduction polarity is controlled by the buoyancy gradient and the strengths of the plates. With no inclined weak zones, imposed orthogonal compression results in the nucleation of a new fault in the weakest plate leading to the young and positively buoyant plate subducting. However, with an inclined weak zone, the buoyancy contrast controls subduction polarity with the most negatively buoyant plate subducting and a self-sustaining subduction regime obtained after ~300 km of imposed shortening.

Interestingly, this situation is also obtained when including an inverted triangular weak zone on top of the transform fault associated with the serpentinization of the crust and mantle.

In non-cylindrical experiments, taking into account the change along strike of plate strength and buoyancy, the capacity of the transform fault to generate a self-sustaining subduction regime is greatly reduced. Subduction initiates simultaneously with opposite polarity at the two extremities of the transform segment and, at depth, a lithospheric tear is produced that separates the two subducting slabs. In the center of the transform fault, the lack of buoyancy or strength contrast between the two plates leads to multiple thrusts with variable polarities, overlapping each other, and each accommodating too little shortening to become the new plate boundary. This indicates that additional mechanical work is required in the center of the transform fault which prevents the establishment of a self-sustaining subduction regime.

The paper is open access and can be obtained here.


Figure 1. Maps of seafloor age (Müller et al., 1997, 2008), plate thickness and buoyancy relative to underlying mantle for fast spreading ridge separating the Pacific and Antarctic plates (A–C), or slow spreading ridge separating the Africa and Antarctic plates (C,E,F). Half spreading rates from GSRM (Kreemer et al., 2003, 2014). Points labeled 1–6 refer to profiles in Figure 2. Plate thickness and relative buoyancy are calculated using the half-space cooling model and the simple plate structure discussed in the text.


Figure 8. Successive stages of Exp. 21. Panels (A–E) show the surface view of the deformed model with PIV vectors, and convergence-parallel horizontal shortening rate. Panel (F) shows the evolution of the force measured at the trailing edge of the plate on the left-hand side. The variation of plate thickness and buoyancy along the strike of the transform fault resulted in generation of multiple faults.

Another way to look at strain

Here is another way to look at strain from PIV analysis of analogue models. After computing all the spatial derivatives of the velocity or incremental displacements field, the 2D small strain tensor can be assembled.

Below I calculated the principal strain direction from the small strain tensor. This is the direction of one of the principal strain, measured from the x-axis.


And then I calculated the values of the maximum principal strain and minimum principal strain. Two horizontal bands appear in the plot. The white cells in the top right and bottom left corners show that this experiment is a sinistral horizontal shear.

The data presented here is not interpolated. Each cell is a data point obtained from image cross-correlation. The maps could be interpolated and be made smoother.



and then we can claculate the maximum shear. We can see where the shear is occuring, where the difference between the maximum and minimum principal strains is largest.


Finally, since we have the directions and values of the principal small strains, we can plot them as crosses. The length of the line is linearly proportional to the magnitude of the principal strain, the orientation of the line is the orientation of the principal strain, and red indicates contraction while blue indicates extension.


A zoom-in shows the directions of S1 and S2 around the developping plastic shear zone fitting what we expect for a sinistral horizontal shear.


I will integrate this into TecPIV rapidly.

Controls on sill and dyke-sill hybrid geometry and propagation in the crust: The role of fracture toughness


Analogue experiments using gelatine were carried out to investigate the role of the mechanical properties of rock layers and their bonded interfaces on the formation and propagation of magma-filled fractures in the crust.

Water was injected at controlled flux through the base of a clear-Perspex tank into superposed and variably bonded layers of solidified gelatine. Experimental dykes and sills were formed, as well as dyke-sill hybrid structures where the ascending dyke crosses the interface between layers but also intrudes it to form a sill.

Stress evolution in the gelatine was visualised using polarised light as the intrusions grew, and its evolving strain was measured using digital image correlation (DIC).
During the formation of dyke-sill hybrids there are notable decreases in stress and strain near the dyke as sills form, which is attributed to a pressure decrease within the intrusive network. Additional fluid is extracted from the open dykes to help grow the sills, causing the dyke protrusion in the overlying layer to be almost completely drained.

Scaling laws and the geometry of the propagating sill suggest sill growth into the interface was toughness-dominated rather than viscosity-dominated. We define KIc* as the fracture toughness of the interface between layers relative to the lower gelatine layer (KIcInt / KIcG). Our results show that KIc* influences the type of intrusion formed (dyke, sill or hybrid), and the magnitude of KIcInt impacted the growth rate of the sills. KIcInt was determined during setup of the experiment by controlling the temperature of the upper layer Tm when it was poured into place, with Tm < 24 °C resulting in an interface with relatively low fracture toughness that is favourable for sill or dyke-sill hybrid formation. The experiments help to explain the dominance of dykes and sills in the rock record, compared to intermediate hybrid structures.

Tectonophysics 698 (2017) 109–120


Dyke-sill hybrid formation, with fluorescent particles in the gelatine illuminated by a thin vertical laser sheet. The intrusion is viewed perpendicular to the dyke strike direction


Photo of dyke-sill hybrid formation. The intrusion is viewed with polarised light, approximately perpendicular to the strike direction of the dyke. Interference colours indicate the evolving distribution and intensity of stress within the gelatine host.


Dyke-sill hybrid formation. The intrusion is viewed looking down and from the side, onto the interface between the gelatine layers. The position of the interface against the tank wall is indicated by the dashed line.
A) A penny-shaped dyke has propagated through the lower gelatine layer and slightly protruded into the upper layer, with two small sills intruding the horizontal interface where it is intercepted by the dyke margins.
B) The dyke protrusion in the upper layer quickly became arrested as the sills grew.
C) The sills joined together within the interface, continued to grow and then coalesced with one margin of the dyke to create the final dyke-sill hybrid structure.

Slab breakoff: insights from 3D thermo-mechanical analogue modelling experiments




The detachment or breakoff of subducted lithosphere is investigated using scaled three-dimensional thermo-mechanical analogue experiments in which forces are measured and deformation is monitored using high-speed particle imaging velocimetry (PIV). The experiments demonstrate that the convergence rate in a subduction zone determine if and when slab detachment occurs. Slow subduction experiments (with scaled convergence rates ∼1 cm yr −1) have lower Peclet numbers and are characterised by lower tensile strength subducted lithosphere, causing detachment to occur when the downward pull force exerted by a relatively short subducted slab is relatively low. When continental collision is preceded by slow oceanic subduction, the subducted lithosphere therefore need not be very long or extremely negatively buoyant to cause detachment because the subducted oceanic lithosphere is hot and weak. Under such conditions detachment may occur sooner after the onset of continental subduction than previously predicted. In contrast, if a collision is preceded by rapid subduction (∼10 cm yr −1), breakoff will be delayed and occur only when the convergence rate slowed sufficiently to thermally weaken the slab and cause its eventual failure. The analogue experiments further confirm that slab detachment occurs diachronously as it propagates along the plate boundary. Stereoscopic PIV reveals a characteristic strain pattern that accompanies the detachment. Horizontal contraction and subsidence (with scaled values up to 1200 m) in the trench and forearc area preceeds the passage of the detachment, which is followed by horizontal extension and uplift (up to 900 m). High-frequency monitoring captures rapid propagation of the detachment along the plate boundary at rates of up to 100 cm yr −1. However rate is not constant and interaction between the slab and lower mantle or opening of a backarc basin in the upper plate can reduce or stop slab breakoff propagation altogether.



Successive side views of the models in Experiment 1 and 2. Experiment 1 (a-e), the subducting lithosphere is pushed by the piston at the constant velocity of 2.5 × 10 −4 m s −1 (equivalent to ∼10 cm yr −1 in nature). The slab becomes vertical due to the negative buoyancy but does not break. It folds when hitting the rigid plate that models the impenetrable lower mantle. Experiment 2 (f–j), The model is identical to Experiment 1 but it is the upper plate that is pushed instead of the lower plate. The model evolution is similar to Experiment 1 until the slab touches the lower mantle. The slab angle reduces in the late stages (dashed line in panel j).



Successive side views of the models in Experiment 3. The model is identical to that employed in Experiment 1 (Fig. 4), but the imposed rate is one order of magnitude lower, 2.5 × 10 −5 m s −1 (equivalent to ∼1 cm yr −1 in nature). Very slow subduction leads to multiple slab detachments at 2283, 4266 and 6420 s. We note that the repeated detachment caused extension in the trailing edge of the upper plate, and a slab graveyard sitting on top of the rigid upper mantle.



Sketch of propagating slab detachment with distribution of surface deformation and uplift. Horizontal contraction and surface subsidence is generated ahead of the breakoff tip, while horizontal extension and uplift follow.



Maps of earthquakes hypocenters along the Aleutian subduction zone (a), and Java-Sumatra-Andaman subduction zone (b), with profiles showing the Wadati-Benioff zone. Earthquakes hypocenters are represented by circles with diameter proportional to magnitude, and color indicating depth (see profiles for color scale). Hypocenters are from EHB catalogue (Engdahl et al., 1998). White arrows represent the convergence vectors calculated using the MORVEL global kinematic model (DeMets et al., 2010). V is the convergence rate (in mm yr −1), θ is the obliquity (angle between the normal to the trench and convergence vector), and Vn is the convergence in the direction of the profile (i.e. convergence corrected from obliquity, in mm yr −1). Topography/bathymetry from Smith and Sandwell (1997).




Slow oblique subduction along the northern branch of the Caribbean subduction zone. Map shows the topography of the trench characterized by a deep through between 65 and 67°W. Convergence from MORVEL global kinematic model (DeMets et al., 2010), is only 19 mm/yr at 64°W (white arrow) and the obliquity is approximately 67°, which yields about 7 mm/yr of normal convergence in the subduction west of 64°W. 3 North-South profiles are plotted showing that this area is also characterized by a deep negative free air anomaly (Sandwell et al., 2014), the peak of which is located the forearc. Based on our experimental results we propose that both the topography and gravity anomalies are caused by an excess downward pull in the subducted lithosphere due to ongoing slab detachment. Topography/bathymetry from Smith and Sandwell (1997), gravity from Sandwell et al. (2014).


Particle Imaging Velocimetry

I have spent some time last year creating a Particle Imaging Velocimetry system to be employed to monitor analogue modelling experiments of tectonics. Basically an image correlation technique is used together with a calibration function to calculate physical displacements and deformation in the models.

My PIV software has been published and is available now. I will post resources soon to help getting started. But first here are some videos of the software. The software has successfully been used last year by my students David Beckett and Maxime Henriquet. It allowed them to quickly become able to quantify their experiments. The example of exported result in the video below is from Beckett’s work. The shear zone setup in the first video is from Maxime’s work.


The material or views expressed on this Blog are those of the author and do not represent those of the University.  Please report any offensive or improper use of this Blog to RPS@newcastle.edu.au.
Skip to toolbar