2-D ﬁnite displacements and strain from Particle Imaging Velocimetry (PIV) analysis of tectonic analogue models with TecPIV

. Image correlation techniques have provided new ways to analyze 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 ﬁnite displacements and ﬁnite strain tensor. We illustrate with synthetic images how the algorithm produces maps of the velocity gradients, small-strain tensor components, incremental or instantaneous principal 5 strains and maximum shear. The incremental displacements can then be summed up with Eulerian or Lagrangian summation, and the components of the 2-D ﬁnite strain tensor can be calculated together with the ﬁnite principal strain and maximum ﬁnite shear. We benchmark the measures of ﬁnite displacements using speciﬁc synthetic tests for each summation mode. The deformation gradient tensor is calculated from the deformed state and decomposed into the ﬁnite rigid-body rotation and left or right ﬁnite-stretch tensors, allowing the deformation ellipsoids to be drawn. The ﬁnite strain has long been the only quantiﬁed 10 measure of strain in analogue models. The presented software package allows producing these ﬁnite strain measures while also accessing incremental measures of strain. The more complete characterization of the deformation of tectonic analogue models will facilitate the comparison with numerical simulations and geological data and help produce conceptual mechanical models.


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, andreference therein).
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 reference therein).Ideally, the evolution of the scaled model is similar to that of the original, but on a much more convenient time scale (e.g.Buckingham, 1914;Hubbert, 1937;Ramberg, 1967;Shemenda, 1994;Merle,  The peak of the correlation map is precisely (i.e., sub-pixel) located using a Gaussian fit.

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.
The flow field was then obtained through an image correlation technique (e.g.Adrian, 1991Adrian, , 2005;;Raffel et al., 2007).In 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 −1 ) (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 (DSLR) to capture sufficient light during the long exposure time.We acknowledge that such 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., 2015Kavanagh et al., , 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 (i.e., 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.

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, 1991Adrian, , 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.

Multipass 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 analyzed using a relatively large interrogation x position (pixel) x position (pixel) (d) multipass 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.
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).

Multipass 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 32 pixels edge length yield spurious variations in the magnitude of the velocity vectors where rigid-body translation is expected (Fig. 4b-c).
A multipass 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-passes approach with a first window of 64 pixels, a second window or 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-passes approach with interrogation windows of 128, 64, 32 and 16 pixels does further reduce 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 multipass 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 multipass approach with an interrogation window decreasing from 64 to 32 (2 passes, Fig. 4e), or from 128 to 32 (3 passes) provided errors of 0.0632 and 0.0622 pixels respectively.We therefore confirm that the multipass approach with window deformation does provide a more accurate capture of the velocity distribution than a single pass with large overlap (Scarano, 2002).
In order to capture the strain distribution within a shear zone, it is necessary to have a vector grid spacing that is at least 1/8th 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 multipass approach with window deformation does indeed facilitate monitoring narrow shear zones since the small window is employed last and informed by the the previous passes.We therefore recommend a strategy that employs at least 2 or 3 passes, 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.

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: 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 extend 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.

2D 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: and The tensor S is the 2D small strain rate tensor.S is derived from D : Its components are often noted ėxx , ėxy , and ė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 2D incremental strain tensor can be obtained from the 2D strain-rate tensor by multiplying it by the time interval.Figure 5 shows ėxx and ėxy , components of the 2D 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.

2D rotation tensor
The anti-symmetric part of the velocity gradient tensor is : A can be written as: The anti-symmetric tensor can be characterised by one 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).

2D 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: while the second relates to the deviatoric strain or shape change: 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.

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 − 1 and e min = S min − 1.
The orientations (θ p and θ p + π/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).The maximum amount of shear is calculated from the principal strains: 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 figure 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 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).

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 a 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 deforma- tion 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(Boutelier et al., , 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 section 3.2.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 10 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 is: 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.

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 5 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 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.

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 maybe 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 noncollinear 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 is: 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).

Benchmark of Lagrangian sum
To benchmark the Lagrangian sum, we created a series of synthetic images were 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-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 using the Lagrangian mode (Fig. 8d-f).Figure 9a shows the profiles of v across the shear zone at successive steps, while figure 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 consist 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.

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: 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.

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 2D the Lagrangian finite strain tensor is : It can be seen that E is symmetrical and contains second-order terms in contrast to the small strain tensor (Eq.5).

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: 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): 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 case the final stage is identical (Fig. 11).The right-stretch tensor U can be obtained from F: and then the rotation tensor can be obtained from F and U: R = FU −1 (29) In 2-D the rigid body rotation can be characterised by a single rotation angle: The left-stretch tensor V can also be obtained from F: 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 1 stretch lines (Ramsay and Huber, 1987).

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 of 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 many particles may be tracked thanks to PIV informing the final PTV step.

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 multipass 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 single pass with larger overlap (Boutelier, 2016).

Figure 1 .
Figure 1.Principles of PIV.(a): Pre-processed 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.

Figure 2 .Figure 3 .
Figure 2. Multipass 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 4 .
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;

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

Figure 6 .
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 emax and emin 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 .
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 5 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 5 frames; (d): Integral of the shear profile versus time measured from PIV (red crosses) and solution (blue line) every 5 frames; (e): Error of the peak cumulative shear strain; (f): Error of the integrated cumulative shear strain.

Figure 8 .
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 4 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 4 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 .Figure 10 .
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 4 frames; (d): Integral of the shear profile versus time measured from PIV (red crosses) and solution (blue line) every 4 frames; (e): Error on the peak cumulative shear strain; (f): Error on the integrated cumulative shear strain.

Figure 11 .
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 are plotted as ellipses.The colour indicates the amount of finite rigid-body rotation obtained via tensor R (in degrees counterclockwise).Positions of the plotted areas are indicated in Fig. 8.

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 I 1e : first invariant of small strain tensor I 2e : second invariant of small strain tensor θ p : Orientation of principal strain e min : Minimum principal strain e max : Minimum principal strain γ max : Maximum shear U E : Component of Eulerian finite displacement (x) V E : Component of Eulerian finite displacement (y) U L : Component of Lagrangian finite displacement (x) V L : 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 Acknowledgements.The authors acknowledge funding for this study from the Australian Research Council via Discovery Project DP14001200.