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

Image correlation techniques have provided new ways to analyse the distribution of deformation in analogue models of tectonics in space and time. Here, we demonstrate, using a new version of our software package (TecPIV), how the correlation of successive time-lapse images of a deforming model allows not only to evaluate the components of the strain-rate tensor at any time in the model but also to calculate the finite displacements and finite strain tensor. We illustrate with synthetic images how the algorithm produces maps of the velocity gradients, small-strain tensor components, incremental or instantaneous principal strains and maximum shear. The incremental displacements can then be summed up with Eulerian or Lagrangian summation, and the components of the 2-D finite strain tensor can be calculated together with the finite principal strain and maximum finite shear. We benchmark the measures of finite displacements using specific synthetic tests for each summation mode. The deformation gradient tensor is calculated from the deformed state and decomposed into the finite rigid-body rotation and left or right finite-stretch tensors, allowing the deformation ellipsoids to be drawn. The finite strain has long been the only quantified measure of strain in analogue models. The presented software package allows producing these finite strain measures while also accessing incremental measures of strain. The more complete characterisation of the deformation of tectonic analogue models will facilitate the comparison with numerical simulations and geological data and help produce conceptual mechanical models.

1 Introduction

The concept of physical similarity rests on the idea that multiple physical systems may share the same underpinning physical laws (Sterrett, 2009, 2017a, b, and references therein).

Figure 1 Principles of PIV. (a) Preprocessed image at time t with illuminated particles; (b) images at time t and t+Δt are split into small interrogation windows (here 128 pixels), and the interrogation windows are cross-correlated; (c) correlation yields a map, indicating how well the interrogation windows correlate with each other when shifted along x and y axes of the image space (here interrogation windows were 32 pixels wide); (d) the peak of the correlation map is precisely (i.e. sub-pixel) located using a Gaussian fit.

In Earth sciences, scaled models have been employed for over a century to test hypotheses on the driving mechanisms of tectonic processes derived from, and constrained by, a variety of geological and geophysical data (e.g. Koyi, 1997; Ranalli, 2001; Graveleau et al., 2012, and references therein). Ideally, the evolution of the scaled model is similar to that of the original but on a much more convenient timescale (e.g. Buckingham, 1914; Hubbert, 1937; Ramberg, 1967; Shemenda, 1994; Merle, 2015). The physical process can thus be monitored more precisely in time and space in the scaled analogue model than in nature, which facilitates its evaluation and interpretation.

In tectonic analogue modelling, the monitoring of model kinematics – local velocities and displacements – provided particularly useful insights (e.g. Cruden et al., 2006; Boutelier and Cruden, 2008; Duarte et al., 2013). In the last two decades or so, one specific monitoring method for model velocities was established as the new experimenter's gold standard: particle imaging velocimetry (PIV) (e.g. Adam et al., 2005, 2013; Schrank et al., 2008; Rosenau et al., 2009; Saumur et al., 2015; Boutelier et al., 2014; Boutelier and Cruden, 2017, 2018; Kavanagh et al., 2015, 2017; Zwaan et al., 2018; Dooley et al., 2018; Rudolf et al., 2019). PIV was developed in the fluid dynamics community to track fluid flow (Adrian, 2005, and references therein). Fluids were seeded with passive tracers, the motions of which were observed with time-lapse photographs. The flow field was then obtained through an image correlation technique (e.g. Adrian, 1991, 2005; Raffel et al., 2007). In tectonic laboratory experiments, PIV is used to monitor the displacements/velocities of markers in or on the model to compute strain or strain-rate components, thus revealing where, how and when deformation occurs in the model.

Because commercial PIV systems are primarily designed for the quantification of rapid fluid flow, these systems usually employ expensive, high-sensitivity but relatively low-resolution digital cameras (<10 Mpixels). The high sensitivity is required to capture the rapid motion of the tracers in the fast-moving fluid at high frequency. The low resolution is a consequence of the required high frequency because it remains technically challenging to move or store a large number of large images at high rates. However, many analogue modelling experiments of tectonics are conducted at very low rates (i.e. <1 mm s⁻¹) (e.g. Corti et al., 2003) or are not rate sensitive at all, and therefore the rate can be set as low as desired (e.g. Hoth et al., 2007). Thus, lower sampling rates can be used, allowing the use of less expensive, higher-resolution (up to 50 Mpixels) consumer-grade digital cameras (DSLRs) to capture sufficient light during the long exposure time. We acknowledge that such an approach is not satisfactory for analogue tectonic models that may include rapid displacements, such as the models of seismic rupture (e.g. Rosenau et al., 2009; Corbi et al., 2013; Rudolf et al., 2019), magma mixing (e.g. Saumur et al., 2015) or the formation of dykes and sills (e.g. Kavanagh et al., 2015, 2017).

For the purpose of facilitating the PIV analysis of analogue models of tectonics, we developed TecPIV, a MATLAB software package allowing, within a single interface, to calibrate and correct the views of the model, correlate the successive images and produce image outputs with incremental displacements and their spatial derivatives (Boutelier, 2016). Indeed, a number of recently published studies using PIV in the tectonics community report incremental (i.e. instantaneous) measures of displacement (e.g. Schmatz et al., 2010; van Gent et al., 2010; Marshak et al., 2019). However, these measures are more akin to GPS-derived velocity fields and calculated strain-rate maps (e.g. Allmendinger et al., 2005; Gupta et al., 2015; Chousianitis et al., 2015; Gunawan and Widiyantoro, 2019) than measures of finite deformation recorded in rocks (Hindle et al., 2002; McQuarrie, 2002; Oncken et al., 2006; Arriagada et al., 2008). For the comparison of laboratory models to field data related to long-term tectonic processes, finite-strain measures are more useful. Measures of finite displacements can also be obtained using the particle tracking velocimetry (PTV) technique (e.g. Toeneboehn et al., 2018) where individual tracers are tracked. However, the technique works best when the number of tracers is limited. Therefore, the information density (i.e. number of measures of displacement) is lower than generally obtained in PIV. An alternative approach to obtain finite displacements and strain in models is to perform a PIV analysis and sum the measured incremental displacements. Here, we update the previously published PIV package through the implementation of the Eulerian and Lagrangian sums of displacements and components of the 2-D finite strain tensor. We detail and benchmark the technique using synthetic images.

Figure 2 Multi-pass PIV with window deformation. (a) First pass (blue area) with 50 % overlap generates the grey vectors as well as the blue vector in the centre; (b) second step with smaller interrogation window (red area in panel a). The results of the first pass are interpolated linearly, or using a spline function, onto the vector grid of the second pass. The interpolated vectors (red vectors) are employed to shift and deform the interrogation prior to correlation. The correlation is then employed to correct the interpolated vectors.

Figure 3 Synthetic static shear zone test at frame 2 (a) and 100 (b). Width of the shear zone: W_sz=512 pixels; y component of the velocity vector: v=4 pixels. Red is the region of interest (ROI); green is the deformed mask.

2 Particle imaging velocimetry

2.1 Principles

In fluid dynamics experiments monitored with PIV, a plane within the fluid is illuminated using a laser sheet (e.g. Adrian, 2005; Raffel et al., 2007). The light scattered by particles in the fluid is recorded on two successive images. Each image is then divided into small interrogation areas, and the local displacement vector between the two images is determined for each interrogation area by cross-correlation (e.g. Adrian, 1991, 2005; Raffel et al., 2007). The algorithm calculates the shift, or displacement vector, for which each interrogation window best matches the next image (Fig. 1). Given the time interval between the two images and the image scaling, the local flow velocity vector can be deduced.

2.2 Multi-pass with window deformation

Our software package now allows performing the correlation in multiple successive steps (passes) with decreasing window size and with window deformation. A first pass is performed, where the image pair is analysed using a relatively large interrogation area and overlap of 50 %. The resulting grid of vectors is interpolated to the smaller spatial resolution of the second pass. This interpolated vector field is then employed to pre-shift and pre-deform the interrogation area before the correlation (Fig. 2). The result of the second correlation is thus a correction of the interpolated first pass. A third and fourth pass can be added following the same principle (Scarano, 2002).

2.2.1 Multi-pass performance

To test how well the new window deformation algorithm captures narrow shear zones, we created a series of synthetic images with a simple static shear zone (see synthetic shear zone test in Boutelier, 2016). The synthetic images are 2000×2000 pixels and include 240 000 particles with a diameter of 3±1 pixels, corresponding to an area proportion of ca. 60 %. The right-hand side of the shear zone is moving upward at a constant velocity, while the left-hand side is static (Fig. 3). The velocity distribution within the shear zone follows the error function. Therefore, the strain-rate profile across the shear zone is a bell-shaped curve. The imposed and measured vertical velocities are compared, and the error is quantified.

For a narrow shear zone with a width of 64 pixels, a single pass with an interrogation area width of 128 pixels and overlap of 50 % captures the rigid-body translations very well in the sense that there is very little noise. The magnitude of the velocity vector is very accurate and the measure is spatially very consistent. However, it fails to capture the width of the shear zone because the vectors are clearly too far apart (Fig. 4a): the effective grid spacing is of the same order as the shear zone width.

A single pass with an interrogation window of 64 or 32 pixels length with overlap of 50 % captures the width of the shear zone more accurately, but an edge length of 64 pixels is still too large a window size, while a 32 pixel edge length yields spurious variations in the magnitude of the velocity vectors where rigid-body translation is expected (Fig. 4b–c).

A multi-pass approach with a first pass with an interrogation window size of 64 pixels, overlap of 50 % followed by a second pass with a window size of 32 pixels, overlap 50 % and linear window deformation captures the shear zone with the precision of the small interrogation window while significantly reducing the noise outside the shear zone (Fig. 4e).

A three-pass approach with a first window of 64 pixels, a second window of 32 pixels and finally a third window of 16 pixels yields a strain profile which captures the shear zone very accurately and does not display severe noise outside the shear zone (Fig. 4g). A four-pass approach with interrogation windows of 128, 64, 32 and 16 pixels further reduces the standard deviation of the residuals but the improvement is too small to be noticed, while the processing time is increased by another 20 %.

The same data set was processed with a single-pass approach with an interrogation window of 64 pixels and an overlap of 75 %. This provided a vector grid with a resolution of 16 pixels (i.e. a vector every 16 pixels), similar to the grids produced by either a single pass with an interrogation area of 32 pixels and overlap of 50 %, or any multi-pass approach with a final pass with an interrogation area of 32 pixels and overlap of 50 %. Using the larger overlap (75 %), the error on the velocity was 0.1045 pixels. This improved over the single pass using an interrogation window of 32 pixels and an overlap of 50 %, which gives an error of 0.124 pixels (Fig. 4c). However, the multi-pass approach with an interrogation window decreasing from 64 to 32 (two passes; Fig. 4e), or from 128 to 32 (three passes) provided errors of 0.0632 and 0.0622 pixels, respectively. We therefore confirm that the multi-pass approach with window deformation does provide a more accurate capture of the velocity distribution than a single pass with large overlap (Scarano, 2002).

Figure 4 PIV results for the narrow static shear zone synthetic test. Shear zone width is 64 pixels; y-axis velocity component of the right-hand side is 2 pixels per frame. (a) Single-pass interrogation area width: w=128 pixels; (b) single pass with w=64 pixels; (c) w=32 pixels; (d) multi-pass with w=128, followed by w=64 pixels; (e) w=64, followed by w=32 pixels; (f) w=32, followed by w=16 pixels; (g) w=64, followed by w=32 and w=16 pixels; (h) w=128, followed by w=64, w=32 and w=16 pixels.

In order to capture the strain distribution within a shear zone, it is necessary to have a vector grid spacing that is at least one-eighth of the shear zone width. However, for narrow shear zones, this may require the use of a small interrogation window which can lead to an increased error because a small window has a less unique distribution of intensity values than a large one. We note that the use of the multi-pass approach with window deformation does indeed facilitate monitoring narrow shear zones since the small window is employed last and informed by the previous passes. We therefore recommend a strategy that employs at least two or three passes and finishes with an interrogation area sufficiently small that the grid spacing can capture the strain distribution in the monitored shear zones.

However, this strategy may need to be modified depending on the scale of the image, particle size and distribution density. There is a minimum width of the interrogation area associated with the physical size of the particles and image scale. It is obvious that to monitor narrow shear zone, very small particles (several times smaller than the shear zone width), and high image scale (many pixels per physical unit of distance) are better than large particles and low image scale. The minimum width of the interrogation area may also depend on the particle distribution density. If the density is too low, a small interrogation area may not have a sufficiently unique distribution of intensity values to be correlated, or it can be correlated anywhere (i.e. no correlation peak). If that is the case, the interrogation area should be increased to include more particles and obtain a more unique distribution of intensity values required for correlation.

2.3 Velocity gradient tensor

A common practice to evaluate where deformation (or rigid-body rotation) is produced in a model is to calculate the spatial derivatives of the velocity vectors. The gradients are calculated for each component using the central difference for interior data points. For example, for the u component of the velocity vectors,

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{\Delta }u(i,j)}{\mathrm{\Delta }x}}={\displaystyle \frac{u(i+\mathrm{1},j)-u(i-\mathrm{1},j)}{\mathrm{2}\mathrm{\Delta }x}},\end{array}$$

where i and j are the indices of the velocity grids corresponding to axes x and y, respectively, and Δx is the homogeneous spacing of the grid along the x axis. Values along the edges of the domain are calculated with a single-sided difference. Therefore, the produced grids of velocity gradients have the same dimension as the velocity grids. The velocity gradients are particularly useful for the analysis of models in which deformation is not confined to the imaged plane. Figure 3 shows an example of a velocity gradient plotted with the velocity vectors for a synthetic static shear zone test. Also displayed are the extent of the region of interest (red rectangle in Fig. 3) and the mask advected with the vectors (green outline in Fig. 3) at the beginning (frame 2) and the end (frame 100) of the test.

2.4 2-D strain-rate tensor

The PIV-derived velocity gradient tensor D composed of the spatial derivatives of the velocity can be split into symmetric (S) and anti-symmetric (A) parts:

$$\begin{array}{}\text{(2)}& \mathbf{D}=\left[\begin{array}{cc}\mathrm{\Delta }u/\mathrm{\Delta }x& \mathrm{\Delta }u/\mathrm{\Delta }y\\ \mathrm{\Delta }v/\mathrm{\Delta }x& \mathrm{\Delta }v/\mathrm{\Delta }y\end{array}\right],\end{array}$$

and

$$\begin{array}{}\text{(3)}& \mathbf{D}=\mathbf{S}+\mathbf{A}.\end{array}$$

The tensor S is the 2-D small-strain-rate tensor. S is derived from D:

$$\begin{array}{}\text{(4)}& \mathbf{S}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathbf{D}+{\mathbf{D}}^T\right),\end{array}$$

or

$$\begin{array}{}\text{(5)}& \mathbf{S}=\left[\begin{array}{cc}\mathrm{\Delta }u/\mathrm{\Delta }x& \frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{\Delta }u/\mathrm{\Delta }y+\mathrm{\Delta }v/\mathrm{\Delta }x\right)\\ \frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{\Delta }v/\mathrm{\Delta }x+\mathrm{\Delta }u/\mathrm{\Delta }y\right)& \mathrm{\Delta }v/\mathrm{\Delta }y\end{array}\right].\end{array}$$

Its components are often noted ${\dot{e}}_{xx}$, ${\dot{e}}_{xy}$ and ${\dot{e}}_{yy}$. Each component of the strain-rate tensor can be plotted as a scalar together with the velocity vectors. We note that since the velocities are obtained by dividing the incremental displacements obtained from the correlation by the time interval between the two correlated images, the 2-D incremental strain tensor can be obtained from the 2-D strain-rate tensor by multiplying it by the time interval. Figure 5 shows ${\dot{e}}_{xx}$ and ${\dot{e}}_{xy}$, components of the 2-D small-strain rate plotted as scalars together with the velocity vectors in the synthetic vortex test. Synthetic images are created where particles are advected according to velocity distribution corresponding to a vortex with a centre at the centre of the image and a radius of 400 pixels. Within the vortex, the magnitude of the velocity vector increases linearly with the distance from the centre. Outside the vortex, the velocity decreases non-linearly with the inverse of the distance to the centre. The positions of the areas with positive and negative values, the symmetries and reached peak magnitudes fit very well the predicted values from the imposed velocity field, although some anomalies are produced along the left and top edges of the image space.

2.5 2-D rotation tensor

The anti-symmetric part of the velocity gradient tensor is

$$\begin{array}{}\text{(6)}& \mathbf{A}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathbf{D}-{\mathbf{D}}^T\right),\end{array}$$

or

$$\begin{array}{}\text{(7)}& \mathbf{A}=\left[\begin{array}{cc}\mathrm{0}& \frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{\Delta }u/\mathrm{\Delta }y-\mathrm{\Delta }v/\mathrm{\Delta }x\right)\\ \frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{\Delta }v/\mathrm{\Delta }x-\mathrm{\Delta }u/\mathrm{\Delta }y\right)& \mathrm{0}\end{array}\right].\end{array}$$

A can be written as

$$\begin{array}{}\text{(8)}& \mathbf{A}=\mathit{\omega }\left[\begin{array}{cc}\mathrm{0}& -\mathrm{1}\\ \mathrm{1}& \mathrm{0}\end{array}\right],\end{array}$$

with

$$\begin{array}{}\text{(9)}& \mathit{\omega }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{\Delta }v}{\mathrm{\Delta }x}}-{\displaystyle \frac{\mathrm{\Delta }u}{\mathrm{\Delta }y}}\right).\end{array}$$

The anti-symmetric tensor can be characterised by a single parameter, the vorticity ω, which can be plotted as a scalar value together with the velocity vectors (Fig. 5c). Our vortex synthetic test shows a constant and homogeneous value of the vorticity within the vortex, as predicted from the imposed velocity distribution (Fig. 5c).

2.6 2-D invariants of deformation

The components of a tensor change under a change of coordinate system associated with a rotation. However, certain tensor properties remain constant regardless of the choice of coordinate system and are thus called invariants (e.g. Jaeger et al., 2007; Allmendinger et al., 2012). In 2-D, there are only two invariants. The first invariant relates to the area strain or size change:

$$\begin{array}{}\text{(10)}& I_{\mathrm{1}e}=\mathrm{tr}\left(\mathbf{S}\right)=e_{xx}+e_{yy},\end{array}$$

while the second relates to the deviatoric strain or shape change:

$$\begin{array}{}\text{(11)}& I_{\mathrm{2}e}=\det \left(\mathbf{S}\right)=e_{xx}{e}_{yy}-e_{xy}^{\mathrm{2}}.\end{array}$$

Our software package now allows plotting these scalar invariants together with the vectors (Fig. 5d). Our synthetic vortex test shows that plotting the second invariant of the small-strain tensor allows highlighting the area where the greatest changes in velocity magnitude are, regardless of the orientation of the velocity vector (Fig. 5d). It is an objective measure of deviatoric strain in the model.

Figure 5 PIV results of the vortex synthetic test illustrating some of the scalar parameters that can be derived from the velocity field: (a) e_xx; (b) e_xy; (c) vorticity; (d) second invariant of the small-strain tensor I_2e.

2.7 Incremental principal strain

The principal strains are two mutually perpendicular directions along which there are no shear strains, only stretch. The stretches in the directions of the principal strains are principal stretches S_max and S_min, and the extensions in these directions are principal extensions $e_{\max }=S_{\max }-\mathrm{1}$ and $e_{\min }=S_{\min }-\mathrm{1}$.

The orientations (θ_p and ${\mathit{\theta }}_{\mathrm{p}}+\mathit{\pi }/\mathrm{2}$) and magnitudes of the incremental principal strains, stretches or extensions can be obtained using the transformation equations describing the changes in the small-strain tensor associated with a rotation of the reference frame (e.g. Jaeger et al., 2007; Allmendinger et al., 2012).

$$\begin{array}{}\text{(12)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{p}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\tan }^{-\mathrm{1}}\left({\displaystyle \frac{\mathrm{2}{e}_{xy}}{e_{xx}-e_{yy}}}\right)\text{(13)}& {\displaystyle }{e}_{\max }& {\displaystyle }={\displaystyle \frac{e_{xx}+e_{yy}}{\mathrm{2}}}+\sqrt{{\left({\displaystyle \frac{e_{xx}-e_{yy}}{\mathrm{2}}}\right)}^{\mathrm{2}}+e_{xy}^{\mathrm{2}}}& {\displaystyle }\text{(14)}& {\displaystyle }{e}_{\min }& {\displaystyle }={\displaystyle \frac{e_{xx}+e_{yy}}{\mathrm{2}}}-\sqrt{{\left({\displaystyle \frac{e_{xx}-e_{yy}}{\mathrm{2}}}\right)}^{\mathrm{2}}+e_{xy}^{\mathrm{2}}}\end{array}$$

The maximum amount of shear is calculated from the principal strains:

$$\begin{array}{}\text{(15)}& {\mathit{\gamma }}_{\max }=e_{\max }-e_{\min }.\end{array}$$

The incremental maximum shear always occurs in a coordinate system that is rotated 45^∘ from the incremental principal coordinate system.

Once the magnitudes and orientation of the principal strains are calculated, the incremental or instantaneous deformation ellipse can be plotted for each point of our PIV vector grids. In Fig. 6, the velocity vectors of a synthetic shear zone test are plotted, and the maximum shear is plotted as a scalar (Fig. 6a). All velocity vectors possess the same component (u) along the x axis, but the component v along the y axis changes across a shear zone parallel to the y axis and located at x=580 pixels at frame 25. Consequently, the measure of γ_max plotted as a scalar value shows the position of the shear zone and distribution of shear within the shear zone (Fig. 6a). Then the principal strains are calculated, and the instantaneous principal extensions are plotted as a divergent (negative extension) or convergent (positive extension) couple of arrows in a subset of the region of interest. A deformation ellipsoid could have been drawn as well, but the incremental deformation between successive images is very small, and the incremental deformation ellipsoids would be very close to circular. It can be seen, however, that a unit circle in the shear zone becomes an ellipse with a minor axis corresponding to the converging arrows and long axis corresponding to the diverging arrows. The orientation of the ellipse associated with the incremental extension does fit the imposed movement (see velocity vectors in Fig. 6a).

Figure 6 Incremental extensions from PIV. (a) Synthetic PIV test with vertical shear zone. Horizontal velocity is homogeneous, but vertical velocity changes sign across a vertical shear at x=580 pixels at frame 25. The shear is illustrated with γ_max, the difference between the principal strains. (b) Principal extensions e_max and e_min across the shear zone. Deformation ellipsoids could have been drawn as well, but the incremental deformation between successive images is very small and the incremental deformation ellipsoids would be very close to circular.

Figure 7 Benchmark of the Eulerian sum of displacements using a synthetic static shear zone. (a) Cumulative displacement profiles across the static shear zone (profile plotted every five frames only for clarity); (b) gradients of the cumulative displacement profiles from panel (a); (c) evolution of the peak of the gradient of the cumulative displacements (i.e. maximum shear, γ) versus time (assuming 1 s between frames in synthetic tests) measured from PIV (red crosses) versus the solution (blue line) every five frames; (d) integral of the shear profile versus time measured from PIV (red crosses) and solution (blue line) every five frames; (e) error of the peak cumulative shear strain; (f) error of the integrated cumulative shear strain.

Figure 8 PIV results for the moving shear zone synthetic test. The u component of the imposed velocity is homogeneous. The v component changes sign across a shear zone that is advected with u. (a, b, c) Velocity vectors with the velocity gradient Δv∕Δx at frames 10, 30 and 60. The location of the shear zone moves to the right over time. Vectors are plotted every four points in each direction for clarity. (c, d, e) Lagrangian sum of displacements at the same stages plotted as deformed grids. The deformed grid is plotted every four points in each direction for clarity. The scalar parameter is the gradient along the x axis of the cumulative displacement along the y axis: ΔV∕Δx. The magnitude of the gradient increases over time, but the shear zone width remains. The red box shows the area plotted in Fig. 12.

Figure 9 Benchmark of the Lagrangian sum of displacements. (a) Cumulative displacement profiles across the moving shear zone (profile plotted every 10 frames only for clarity); (b) gradients of the cumulative displacement profiles from panel (a); (c) evolution of the peak of the gradient of the cumulative displacements (i.e. max shear, γ) versus time (assuming 1 s between frames in synthetic tests) measured from PIV (red crosses) versus the solution (blue line) every four frames; (d) integral of the shear profile versus time measured from PIV (red crosses) and solution (blue line) every four frames; (e) error on the peak cumulative shear strain; (f) error on the integrated cumulative shear strain.

Figure 10 Illustration of incremental and cumulative Lagrangian displacements and track length. Material point (red dot) moves from position 0 to 5. The velocity vectors at time step i are u_i,v_i. The cumulative displacements are U_i,V_i (dashed line). The track length is the length of the path: the length of the green line.

Figure 11 Polar decomposition of the deformation gradient tensor F into a finite rigid-body rotation tensor R and either a left-stretch tensor V or right-stretch tensor U. I is the identity matrix. Note that F and V have the same orientations of their principal components. The green and orange lines are initially parallel to the axes of the Cartesian reference frame and here illustrate the rotation of the unit circle (top). The dashed lines are the long and short axes of the deformation ellipsoids.

Figure 12 Deformation gradient tensor F and right-stretch tensor U calculated from the Lagrangian sum of displacements at frame 30 and 60 of the moving shear zone synthetic test. For clarity, only one in four points is plotted as an ellipse. The colour indicates the amount of finite rigid-body rotation obtained via tensor R (in degrees anticlockwise). Positions of the plotted areas are indicated in Fig. 8.

3 Finite displacements and strain

The calculation of finite deformation is very useful for the analysis and interpretation of tectonic models because it is finite deformation that is most commonly and easily observed in the field. In order to calculate the finite deformation in the model, one may first calculate the sum of incremental displacements. This sum can be performed using an Eulerian or Lagrangian approach.

The Eulerian description of a flow field focuses on specific locations in space through which the fluid flows (e.g. Batchelor, 1967). The Eulerian sum of displacements is therefore the measure of the cumulative flow through these specific locations. On the other hand, the Lagrangian description of a flow field defines the flow as number of individual fluid parcels that move through space and time. The Lagrangian sum tracks the cumulative displacements of each parcel or particle in the fluid (e.g. Jaeger et al., 2007; Allmendinger et al., 2012).

The Eulerian approach is commonly used to track viscous fluid flow where the mechanical properties are not strain dependent. For more solid-like behaviour, the Lagrangian approach allows calculating the cumulative strain even when the deformation structures (shear zones) are not fixed in space. For example, analogue models of shortened lithosphere commonly contain conjugate shear zones that are translated with progressive shortening (Boutelier et al., 2012, 2014). An Eulerian approach, however, may be sufficient if the assumption that the deformation structures are fixed in space can be reasonably made (for example, when it is controlled by the experimental setup). An Eulerian sum of displacements on a progressively translated shear zone would produce smearing of the cumulative strain over the various spatial positions of the shear zone over time, as elucidated further in Sect. 3.2.1.

3.1 Eulerian sum of displacements

The Eulerian sum is straightforward when the employed region of interest (ROI) and mask are static since the displacements can be summed up at points, which are fixed in space throughout the entire recording. When the mask and ROI deform with the flow, the positions of the points where incremental displacements are known will change, and the new points need to be interpolated. For simplicity, we perform the Eulerian sum only within a fixed mask.

The components U_E,V_E of the finite displacement at coordinates x,y and time t_e using the Eulerian approach are

$$\begin{array}{}\text{(16)}& {\displaystyle }{U}_{\mathrm{E}}(x,y,t_{\mathrm{e}})& {\displaystyle }=\sum _{\mathrm{0}}^{t_{\mathrm{e}}}u(x,y,t)\mathrm{\Delta }t\text{(17)}& {\displaystyle }{V}_{\mathrm{E}}(x,y,t_{\mathrm{e}})& {\displaystyle }=\sum _{\mathrm{0}}^{t_{\mathrm{e}}}v(x,y,t)\mathrm{\Delta }t,\end{array}$$

with $u(x,y,t)$ and $v(x,y,t)$ being the components of the velocity vector obtained by image correlation at coordinates x,y and time t, and Δt is the time interval between correlated images.

3.1.1 Benchmark of Eulerian sum

To benchmark the Eulerian sum, we employed synthetic images with a static shear zone. The images are 1000×1000 pixels and include 60 000 particles with a diameter 3±1 pixels. The u component of the velocity is 0, and the v component varies across a shear zone parallel to the y axis at coordinate x=500 pixels. The shear zone width is 512 pixels, and the velocity varies in the shear zone following the error function from 0 on the left to 4 pixels per frame on the right. Figure 3 shows two stages of the benchmark test at the beginning and end of the simulation. Figure 7a shows profiles of the Eulerian cumulative displacement across the shear zone every five frames during the simulation. Figure 7b shows the gradient of the cumulative displacement along the same profiles (i.e. cumulative shear strain). It can be seen that the shear zone remains static but the amount of accommodated shear increases over time. We compared the peak value of the cumulative shear in the profile to the value derived from the velocity field employed to advect the particle (Fig. 7c). We also compared the integrated cumulative shear in the profile to the values derived from the velocity field employed to advect the particles. Both peak values of the cumulative shear and the integrated cumulative shear, respectively, increase linearly with rates very similar to the expected behaviour. We quantified the difference between the values derived from the PIV analysis and the solution. For both the peak cumulative shear and the integrated cumulative shear, the error is largest at the beginning of the sum and reduces to about 2 % for the peak shear strain and less than 0.5 % for the integrated cumulative shear strain. We conclude that the Eulerian sum provides a very reliable way to estimate the amount of shear accommodated in the static shear zone.

3.2 Lagrangian sum of displacements

The Lagrangian sum requires that the incremental displacements, which are summed up, be advected with the flow. This poses some difficulties because the displacement is defined at each time step only on the points of the vector grid but not in between. However, the incremental displacement (i.e. between successive images) is generally significantly smaller than the vector grid resolution. For example, the vector grid may have points every 32 pixels but the incremental displacement may be up to a few pixels only. The vectors that are summed up are thus vectors which have been interpolated from neighbouring points on the vector grid.

Here, the value at each requested point that falls in between known vector points is calculated using the three nearest non-collinear neighbours. If the three nearest neighbours method fails because the points are collinear, the value of the unsampled point is calculated as the weighted average of the known values of the N nearest neighbours. The weights of the neighbours are inversely related to the distances between the prediction location and the sampled locations. This method allows rapid calculation of the Lagrangian sum of displacements even for large grids.

The components U_L,V_L of the finite displacement of a material point a at time t_e are

$$\begin{array}{}\text{(18)}& {\displaystyle }{U}_{\mathrm{L}}(a,t_{\mathrm{e}})& {\displaystyle }=\sum _{\mathrm{0}}^{t_{\mathrm{e}}}u\left(X\right(a,t\left)\right)\mathrm{\Delta }t\text{(19)}& {\displaystyle }{V}_{\mathrm{L}}(a,t_{\mathrm{e}})& {\displaystyle }=\sum _{\mathrm{0}}^{t_{\mathrm{e}}}v\left(X\right(a,t\left)\right)\mathrm{\Delta }t,\end{array}$$

where X(a,t) is the location of the material point a at time t, and u(X(a,t)) is the interpolated velocity component of the material point a at location X(a,t).

3.2.1 Benchmark of Lagrangian sum

To benchmark the Lagrangian sum, we created a series of synthetic images where the particles are moved homogeneously along the x axis but heterogeneously along the y axis. On the left-hand side of a shear zone parallel to the y axis, the particles move downward, while on the right-hand side of the shear zone the particles move up. In between is a shear zone where the y-axis velocity component varies along x as an error function. The shear profile across the shear zone is thus a bell-shaped curve. At each time step, the x position of the shear zone is advected along the x axis with the particle. An Eulerian sum would thus smear the shear zone over the various positions it has in the recording, but a Lagrangian sum should maintain the width of the shear zone. The images included 240 000 particles with a diameter of 3±1 pixels. The width of the shear zone is 512 pixels, the horizontal incremental displacement is 3 pixels, while the y-axis velocity component varies between −1.5 and 1.5 pixels. We calculated the incremental displacements (Fig. 8a–c) and summed them using the Lagrangian mode (Fig. 8d–f). Figure 9a shows the profiles of v across the shear zone in successive steps, while Fig. 9b shows the gradient of the finite displacement component V. It can be seen that the peak is shifting appropriately, and the width of the shear zone appears to be maintained properly. To quantify further, we calculated the peak of the shear profile and its integral at each time step and compared them with values computed from the equation employed to advect the particles in the synthetic images (Fig. 9c–d). The error on the magnitude of the peak increases over time but is less than 0.1 % after 100 steps, corresponding to a relative displacement along y of 300 pixels. The error on the integrated shear decreases rapidly in the first steps and then continues to slowly decrease until step 100. The final error on the integrated shear strain is less than 1 %.

We conclude that the calculation of the Lagrangian sum is precise when the incremental displacement can be measured accurately. Of course, if the quality of the experimental images does not permit a precise measure of the incremental displacements, the finite displacements will also reflect this inaccuracy. However, the summation process is itself accurate. Therefore, if the inaccuracy in the measure of the incremental displacements consists of random error on the displacements, the cumulative sum of displacements might increase the strength of the signal. However, if the error is systematic (i.e. same in every frame) due to, for example, a poor or wrong image calibration, then the cumulative sum may amplify this condition and be unusable. We note here that we employed a relatively simple linear interpolation from the three nearest neighbours to calculate the interpolated velocity and incremental displacement that are summed up. This interpolation technique was chosen for its rapidity and ability to get values at precisely the location of the unsampled point. It is possible that a spline or cubic polynomial employing more known neighbouring points might provide an even better accuracy.

3.2.2 Track length

The finite strain tensor is calculated using the Lagrangian cumulative displacement above and therefore compares the initial and final states, irrespective of the strain path (Fig. 10). A measure of the length of this path is provided with the following parameter:

$$\begin{array}{}\text{(20)}& L(a,t_{\mathrm{e}})=\sum _{\mathrm{0}}^{t_{\mathrm{e}}}\mathrm{\Delta }t\sqrt{u(a,t{)}^{\mathrm{2}}+v(a,t{)}^{\mathrm{2}}}.\end{array}$$

The components of the incremental displacement of a material point a are interpolated and employed to calculate the magnitude of the incremental displacement, which is summed over the duration of the recording.

3.3 Finite strain tensor

The Lagrangian finite strain tensor E is calculated using the gradients of the Lagrangian cumulative displacements and the initial relative position of the points. In 2-D, the Lagrangian finite strain tensor is

$$\begin{array}{}\text{(21)}& \mathbf{E}=\left[\begin{array}{cc}{E}_{xx}& E_{xy}\\ E_{xy}& E_{yy}\end{array}\right],\end{array}$$

with

$$\begin{array}{}\text{(22)}& {\displaystyle }{E}_{xx}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }x}}+{\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }x}}+\left({\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }x}}{\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }x}}+{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }x}}{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }x}}\right)\right)\text{(23)}& {\displaystyle }{E}_{yy}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }y}}+{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }y}}+\left({\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }y}}{\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }y}}+{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }y}}{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }y}}\right)\right)\text{(24)}& {\displaystyle }{E}_{xy}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }y}}+{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }x}}+\left({\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }x}}{\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }y}}+{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }x}}{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }y}}\right)\right)\text{(25)}& {\displaystyle }{E}_{yx}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }x}}+{\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }y}}+\left({\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }y}}{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }x}}+{\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }y}}{\displaystyle \frac{\mathrm{\Delta }U}{\mathrm{\Delta }x}}\right)\right).\end{array}$$

It can be seen that E is symmetrical and contains second-order terms in contrast to the small-strain tensor (Eq. 5).

3.4 Rotation and stretch tensors

With finite strain, the displacement gradient tensor cannot be decomposed additively into symmetric strain and anti-symmetric rotation parts. However, the deformation gradient tensor can be decomposed into the product of two tensors: one accounting for the stretch and one for the rigid-body rotation (e.g. Allmendinger et al., 2012).

The deformation gradient tensor is defined as the gradient of the deformed distances of material points relative to their initial distances:

$$\begin{array}{}\text{(26)}& F_{ij}={\displaystyle \frac{\mathrm{\Delta }{X}_i}{\mathrm{\Delta }{x}_j}},\end{array}$$

with ΔX_i the distance between the points in the direction i (x or y) in the deformed state, while Δx_j is the distance between the points in the direction j in the undeformed state. Since we can extract the position of points in the deformed state using the Lagrangian sum of displacements and already know the initial positions of these points, the deformation gradient tensor can be extracted from the Lagrangian sum.

The deformation gradient tensor can be decomposed multiplicatively in two ways (e.g. Allmendinger et al., 2012):

$$\begin{array}{}\text{(27)}& \mathbf{FI}=\mathbf{RUI}=\mathbf{VRI},\end{array}$$

where R is the rotation tensor, U is the right-stretch tensor, V is the left-stretch tensor, and I is the identity matrix. In the first case, the rigid-body rotation occurs first, followed by stretching, while in the second case stretching occurs before rigid-body rotation. In both cases, the final stage is identical (Fig. 11). The right-stretch tensor U can be obtained from F:

$$\begin{array}{}\text{(28)}& \mathbf{U}={\left({\mathbf{F}}^T\mathbf{F}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

and then the rotation tensor can be obtained from F and U:

$$\begin{array}{}\text{(29)}& \mathbf{R}={\mathbf{FU}}^{-\mathrm{1}}.\end{array}$$

In 2-D, the rigid-body rotation can be characterised by a single rotation angle:

$$\begin{array}{}\text{(30)}& {\mathrm{\Theta }}_{\mathrm{f}}={\cos }^{-\mathrm{1}}\left(R_{\mathrm{11}}\right).\end{array}$$

The left-stretch tensor V can also be obtained from F:

$$\begin{array}{}\text{(31)}& \mathbf{V}={\left({\mathbf{FF}}^T\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

In our code, the components of F are calculated for each Lagrangian point, and then F is assembled. Next, U, R and V are calculated. Finally, the magnitudes and orientations of the principal values are calculated for F and U.

Figure 12 shows the finite strain ellipses using the deformation gradient tensor F and the right-stretch tensor U at frames 30 and 60 of our synthetic test on a moving shear zone employed to measure the accuracy of the Lagrangian sum (Figs. 8 and 9). As expected, the ellipsoid becomes more elongated near the centre line of the shear zone and with time. The amount of finite rigid-body rotation also increases towards the centre line of the shear zone and over time. The measure of stretch before rigid-body rotation from U shows the principal stretches forming the expected S₁ stretch lines (Ramsay and Huber, 1987).

4 Future developments

Future developments of the software package may include introducing an adaptive refinement of the interrogation window using a measured rigid-body movement. Where incremental displacements correspond to a rigid body, no increase in resolution is required, while where a “significant” departure from rigid-body displacements is observed, a refinement is required to capture the strain distribution. This would be similar to getting the peak of Fig. 4 with high resolution (e.g. Fig. 4f–h) but the base with very robust but sparse low resolution (Fig. 4a,b). This idea, however, requires splitting the ROI into a number of ROIs with different resolutions.

Another future development is blending PIV and PTV techniques to get a super-resolution PIV (e.g. Stitou and Riethmuller, 2001; Ha et al., 2012; Korevaar et al., 2015). PTV usually works best when the number of particles is small. This is generally the case within the interrogation window after several refinements. Eventually, the interrogation window becomes so small that it does not contain enough particles to perform the PIV correlation. However, the software could then switch to a PTV mode and track the few particles within the interrogation windows. This would produce the highest possible resolution for flow monitoring as well as high vector density since all the particles may be tracked thanks to PIV informing the final PTV step.

5 Conclusions

Particle imaging velocimetry of analogue models of tectonics can be performed with inexpensive, high-resolution commercial DSLR cameras and open-sourced scripts for MATLAB. Our GUI package TecPIV provides a way to perform the calibration and rectification of the model images prior to the correlation of successive image pairs.

A new correlation algorithm has been implemented, allowing the use of the multi-pass procedure with decreasing window size and window deformation. Synthetic tests confirm that this method improves the ability to capture displacements in and around narrow deformation bands compared to the previous approach using a single pass with larger overlap (Boutelier, 2016).

A region of interest is defined as a rectangular box encompassing the deforming model, inside which vectors are calculated. A deforming mask is also calculated for each time step. The mask corresponds to the ROI for the first time step but is advected with the vectors in subsequent steps.

The software package now provides access to the 2-D small-strain tensor components, its invariants as well as the principal extensions, and maximum shear. The incremental deformation ellipsoid, using principal stretches instead of extension, can be derived everywhere in the imaged model and for each time step.

Two approaches are proposed for the calculation of the finite displacements: Eulerian or Lagrangian sums. The Lagrangian sum may be more appropriate for the analysis of tectonic models, as it allows calculating the amount of deformation accommodated by a moving feature. Synthetic tests demonstrate that when the incremental displacements are accurate, the sums are also accurate, as only a very small error is introduced due to the interpolation required for the summation process. The Lagrangian sum of displacements allows the calculation of the 2-D Lagrangian finite strain tensor and its invariants. In addition, the 2-D deformation gradient tensor is calculated from the model deformed state (i.e. from the Lagrangian sum of displacements and initial positions). The 2-D deformation gradient tensor can then be decomposed into the finite rigid-body rotation tensor and either a left or right finite-stretch tensor. Deformation ellipses can be plotted and the effect of the rigid-body rotations and stretch measured and visualised.

We conclude that the PIV analysis of 2-D plane-strain analogue models of tectonics can provide simultaneously the incremental and cumulative strain tensor. The finite deformation was traditionally analysed in analogue models using passive markers (i.e. Graveleau et al., 2012, and references therein). The particle imaging velocimetry technique born of fluid dynamics experiments has recently provided access to incremental, or instantaneous, deformation rates but often neglected finite strain. The cumulative effects of finite strain are required to permit better comparisons of models to geological field observations. The calculation of the Lagrangian sum allows bringing back the measure of the finite deformation.

Appendix A: Nomenclature

x	Coordinate along x axis
y	Coordinate along y axis
u	Component of velocity vector along x axis
v	Component of velocity vector along y axis
v	Velocity vector
s	Displacement vector
Δt	Time between images
Δx	Resolution of vector grids along x axis
Δy	Resolution of vector grids along y axis
i, j	Indices of grids along x and y dimensions
D	Velocity gradient tensor
S	Symmetric part of velocity gradient tensor
A	Anti-symmetric part of velocity gradient tensor
ω	Vorticity
I1e	First invariant of small-strain tensor
I2e	Second invariant of small-strain tensor
θp	Orientation of principal strain
emin	Minimum principal strain
emax	Minimum principal strain
γmax	Maximum shear
UE	Component of Eulerian finite displacement (x)
VE	Component of Eulerian finite displacement (y)
UL	Component of Lagrangian finite displacement (x)
VL	Component of Lagrangian finite displacement (y)
L	Lagrangian track length
E	Lagrangian finite strain tensor
F	Deformation gradient tensor
R	Finite rigid-body rotation tensor
U	Right-stretch tensor
V	Left-stretch tensor
Θf	Finite rigid-body rotation angle