Mapping and classifying large deformation from digital imagery: application to analogue models of lithosphere deformation

Particle Image Velocimetry (PIV), an image cross-correlation technique, is widely used for obtaining velocity fields from series of images of deforming objects. Rather than instantaneous velocities, we are interested in reconstructing cumulative deformation, and use PIV-derived incremental displacements for this purpose. Our focus is on analogue models of tectonic processes, which can accumulate large deformation. Importantly, PIV provides incremental displacements during analogue model evolution in a spatial reference (Eulerian) frame, without the need for explicit markers in a model. We integrate the displacements in a material reference (Lagrangian) frame, such that displacements can be integrated to track the spatial accumulative deformation field as a function of time. To describe cumulative, finite deformation, various strain tensors have been developed, and we discuss what strain measure best describes large shape changes, as standard infinitesimal strain tensors no longer apply for large deformation. PIV or comparable techniques have become a common method to determine strain in analogue models. However, the qualitative interpretation of observed strain has remained problematic for complex settings. Hence, PIV-derived displacements have not been fully exploited before, as methods to qualitatively characterize cumulative, large strain have been lacking. Notably, in tectonic settings, different types of deformation extension, shortening, strike-slip can be superimposed. We demonstrate that when shape changes are de-


Introduction
Determination of the cumulative strain of an object provides insight in the total shape change resulting from prolonged deformation, and aids in the characterization of the type of deformation leading to deformed structures. A technique such as Particle Image Velocimetry (PIV) allows, by cross-correlating subsequent images, to obtain the incremental displacement field in between two time steps (e.g. Thielicke and Stamhuis [2014]). From these displacements, incremental deformation can be derived. In a solid earth context, image cross-correlation can be used to observe displacements or deformation in physical analogue models of longterm tectonics. Applications range from mapping surface motions [Hampel et al., 2004], to the determination of strain in the horizontal plane [Boutelier and Oncken, 2011] or in side view [Cruz et al., 2008], to the imaging of mantle flow in models [Funiciello et al., 2006]. Similar correlation techniques are employed to detect near step-wise displacements from coseismic deformation using optical satellite imagery [Kääb et al., 2017, Sotiris et al., 2018, SAR data [Morishita et al., 2017] or a combination of both [Lauer et al., 2020].
As PIV provides displacements in an Eulerian (i.e. space-fixed) reference frame, to integrate large cumulative deformations we need to convert these displacements to a Lagrangian (material-fixed) frame. Namely, we are primarily interested in the deformation of material, rather than the deformation at a location in space through which material advects. One can use the Eulerian displacements to follow material, provided that PIV-derived displacements are continuously available in time. Senatore et al. [2013] and [Stanier et al., 2016] show that PIV analysis of deformation of granular material allows for high resolution mapping of cumulative strain. Recently, Boutelier et al. [2019] illustrated the potential of integrated PIV displacements based on synthetic PIV results for simple shear zones, using appropriate finite strain descriptions.
As an alternative to PIV, displacement detection methods have been used that rely on tracing passive markers, manually [Schellart et al., 2003], or digitally [Fischer and Keating, 2005, Cardozo and Allmendinger, 2009, Duarte et al., 2013, Guillaume et al., 2013, Chen et al., 2015; for a review see section 5.3 of Schellart and Strak [2016]. Tracking of passive markers provides Lagrangian displacements directly, but faults can be difficult to observe due to limited spatial resolution [Nilforoushan et al., 2008]. PIV relies on the presence of contrasting texture of the study object in view, and allows for a higher spatial resolution of the deformation field compared to methods were individual markers are traced through time. The final resolution depends on the contrast of the surface texture and image resolution [White et al., 2003]. PIV or similar auto-correlation techniques (digital volume correlation) can also be used to map displacements in 3D, for example using X-ray tomography to visualize internal deformation of analogue models [Adam et al., 2013, Poppe et al., 2019 or experimentally deformed rocks [Lenoir et al., 2007, Mao andChiang, 2016]. For transparent media, stereoscopic PIV has also proven itself to image 3D velocity fields in subduction models [Strak and Schellart, 2014]. However, in our study we focus on 2D deformation occurring in a horizontal plane, as our 3 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International study materials are not transparant and we use planar view images to determine deformation.
In this work we seek to infer how the different types of tectonic deformation -extension, shortening and strike-slip -contribute to the cumulative deformation, documented through physical analogue modeling. A PIV-derived classification of the type of strain can so complement visual inspections of the type of faulting (e.g., as in Willingshofer and Sokoutis [2009]), as well as determine distributed deformation in the absence of clear discrete faults.
When strains are small and rotations are negligible, such as holds for interseismic deformation measured using geodetic observations, the type of strain can be obtained using the ratio between the principal strains, such as used by Kreemer et al. [2014] for the interpretation of current day global deformation rates. For large strains this ratio no longer applies, as we will show later. Studies have used derived quantities as maximum shear strains [Cruz et al., 2008] or vorticity [Adam et al., 2005] to resolve shear strain in vertical model crosssections, or the first strain invariant to determine structures that deform under extension or shortening [Galland et al., 2016]. To our knowledge, no existing method allows the simultaneous discrimination of extension, shortening and strike-slip (simple shear), or combinations of those.
In this paper we will discuss various deformation tensors that are common, infinitesimal strain, or less common in geophysical modeling literature, finite strain and finite stretch. We will argue which tensors are best suited to describe large deformation and classify deformation. Based on a finite stretch tensor, we formulate a classification of deformation that can properly discern distributed or localized strike-slip, extension or shortening, or oblique combinations: transpression (strike-slip and orthogonal shortening) or transtension (strikeslip and orthogonal extension). Using both synthetic and analogue models of tectonic deformation, we will show that it is possible to create high resolution maps of finite strain and the type of strain.

Deformation analysis 2.1 Material displacements from time variable spatial displacements
Given incremental displacement vector fields in time, in our case obtained by particle image velocimetry (PIV) [Thielicke and Stamhuis, 2014], we want to determine the deformation field as a function of time. PIV provides the displacements in an Eulerian, space-fixed, undeformed, reference frame, which we want to express in a Lagrangian, material-fixed, deformed reference frame. We define initial material points, from which we reconstruct particle paths through time, using the Eulerian displacements. The left panel of figure 1 shows these initial material points. For the following time steps, material has advected and these material points will no longer coincide with the PIV-derived Eulerian displacement points, see the middle panel. Similar as Senatore et al. [2013], Stanier et al. [2016], Boutelier et al. 4 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Figure 1: Interpolation of Eulerian displacements to follow material and construct a Lagrangian displacement field. Left panel: PIV-derived Eulerian displacements (black arrows) and interpolated displacements (green arrows) on an initial material grid (black) at the first epoch. Middle panel: PIV-derived Eulerian displacements (black arrows) and initial grid (black) and interpolated displacements on deformed material grid (green) at second epoch. Right panel: particle paths after a number of epochs, deforming, rotating and translating the initial grid (black) to the current deformed grid (green).
The current location x of a point at time t that initially had the reference position X we construct as: where δu(X, t) is the incremental Eulerian displacement at time t of the point that occupied position X at the reference epoch. We obtain δu(X, t) by spatial bi-linear interpolation, as δu is provided in a discrete interval by PIV. Figure 1 shows how the original material frame has been deformed by incremental displacements, and how interpolated displacements can be used to construct particle paths. It is clear how for large deformation (i.e., translation, distortion, rotation) the Eulerian displacements quickly deviate from the material displacements. A higher spatial resolution of the grid compared to the PIV displacement field can be useful to mitigate progressive resolution loss in regions with large cumulative extension.
The procedure assumes that the PIV-derived velocities do not contain outliers. In this study we remove outliers based on a threshold displacement value as a function of the standard deviation of all displacement magnitudes [Thielicke and Stamhuis, 2014], followed by a man-5 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International ual check of remaining outliers.

Deformation measures
Throughout this manuscript we use strain and stretch as measures to describe shape changes. As a simple 1D scalar measure, extensional strain is: with initial length L i and final length L f . Extensional strain may stand here for both positive strain, as well as negative, shortening, strain. The range of possible 1D strain values is [-1 ∞], with 0 denoting no deformation. In dimensions larger than 1 multiple definitions are possible for strain, depending on the reference frame and the magnitude of deformation (sections 2.5.1 and 2.5.2). Stretch Λ is a related measure that we will use frequently in this manuscript and in 1D it is defined as: where the stretch is always positive, and a stretch of 1 denotes no deformation. For purely extensional deformation stretch relates to strain as λ = 1 + . Finally, moving to 2D, we will make use of the area ratio dA/dA 0 (i.e., dilatation from a 2D perspective), which is a function of the principal stretches λ 1 , λ 2 as: Finally, engineering shear strain γ is related to the change in angle ψ of two, originally, perpendicular lines: γ = tan(ψ)

Deformation gradient
From the Lagrangian positions in time we derive the deformation gradient F, that we will use later to derive all deformation measures such as strain and stretch. Figure 2 shows the relation between old and new element vectors (i.e., any line within a small unit of material) and displacements. The deformation gradient relates initial vectors of the undeformed material dX to the deformed, final material vector dx, by [Malvern, 1969]: This means that the components of deformation gradient tensor are defined as: 6 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Figure 2: The deformation gradient F relates the deformed state to the undeformed state. A material vector dX that runs from point P to Q (i.e, a side of a small material element) in the original configuration, to the material vector in the deformed configuration dx, with the material that first inhibited points P and Q now in points p and q. Displacement vectors u(P ) and u(Q) provide the displacements between the two states. non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International We can write the displacement u as the difference between current and reference coordinates: Thus, the deformation gradient can be reformulated as a function of the displacements: Here ∂u ∂X or ∇u is the gradient of the displacement vector field, and I is the identity matrix, which represents the current state. This means that the deformation gradient is a simple relation of the partial derivative of the displacement with respect to the reference coordinates. In 2D the deformation gradient becomes: When we regard the deformation as a product of rigid body rotation and translation and distortion, F captures all except the rigid body translation. Namely, rigid body translations do not lead to gradients in displacements.

Computation of the deformation gradient
In equation 1 we discussed the coordinates of points of which we have tracked the change in coordinates x. For the following, we define quadrilateral elements (that in our implementation starts as squares in the undeformed state) and we apply the computation of F to these elements. We make use of linear shape functions to determine the partial derivatives of the displacement ∂u ∂X , see appendix A.

General deformation described with the deformation gradient
Any deformation that can be described by F, can be written as a serial combination of a simple shear, simple extensions along the x or y axes and a rotation [Srinivasa, 2012]: DistortionF contains all shape changes and Q is an orthogonal rotation matrix (Q T Q = I). The distortion can very conveniently be written as a successive multiplication of first a shearing motion followed by a biaxial extension [Freed et al., 2019]: 8 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Figure 3 shows an example of the resulting deformation with applying distortion and rotation. For example, shearing (simple shear) parallel to x can be achieved by engineering strain γ = 0. Subsequently, extension or shortening along x can be achieved with stretch parameter a > 1 and 0 < a < 1, respectively. Similarly extension along y can be achieved using stretch parameter b = 1. Finally, rotation matrix Q rotates the cartesian reference frame x under a counterclockwise angle θ to the new coordinates x : x = Qx. With the rotation given by: We will use the formulation of the distortionF in section 2.7 to classify large shape changes (distortion) as combinations of shearing and extensions.

Strain and stretch tensors
Different measures have been derived to uniquely decompose the deformation described by F into rotations and shape changes. We will discuss a few of these and argue which of these tensors are best equipped to quantitatively describe large changes in shape.

Small deformation: infinitesimal strain
For small deformations, the product of successive deformations (see equations 11 and 12) can be approximated by a sum of successive deformations. In that case, we can write the deformation gradient as the sum of the current state I, shape changes and rotation ω [Allmendinger et al., 2011]: Shape changes are expressed in terms of strain, a quantity closely related to the 1D extension , equation 2. When rotations ω -due to rigid body rotations as well as due to shearare small, we can describe shape changes using the Lagrangian infinitesimal strain tensor (in 2D) [Malvern, 1969]: with normal extension x and y on the diagonals and the shear term xy on the off-diagonals. The infinitesimal strain tensor can be written as a function of the displacement tensor ∇u, and thus can be easily computed from the deformation gradient (F = I + ∇u) as: Complementary to the strain tensor is the rotation tensor ω, which relates to the displacement tensor and strain as ∇u = + ω. For straight material lines in deforming elements, non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International incremental deformation cumulative deformation non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International the rotation differs depending on the orientation and ω represents the average rotation, integrated over all orientations of a small piece of material and is defined as: which can also be written as a function of the deformation gradient:

Example: small strains, rotation and shear
The expression for the strain only makes sense as long as rotations are small, as rotations will affect the diagonal elements of , even in the absence of shape changes. For an object under rigid body rotation with counterclockwise angle θ (a = 1, b = 1, γ = 0), the deformation gradient is: Even in the absence of shape changes the infinitesimal strain then becomes non-zero, which demonstrates the physical insignificance of to describe shape changes when cos(θ) is no longer approximately 1: A second case: for a material subject to simple shear, with shear parallel to x (a = 1, b = 1, γ > 0, θ = 0), the deformation gradient is: and the infinitesimal shear becomes: Material that originally is parallel with the y-axis will elongate, due to the shear parallel to x, see the left panel of figure 3. However, this effect is absent in the infinitesimal strain, as y is zero. For small strains this elongation along y may be small, and can thus be neglected. But it will become significant for large shears, as section 2.5.4 shows. For objects undergoing large deformation, the infinitesimal strain description can thus only be used for small cumulative rotations and small shears. Notably, the summation of incremental strains is only possible when the principal strain axes do not change during deformation [Malvern, 1969], which again rules out application for materials subject to shear or rigid body rotations.
non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International

Large deformation: finite strain
For large deformations there are different strain definitions that are insensitive to rotations, and include extension resulting from shear on diagonal components of strain. Finite strain is related to the squared length change, see for a derivation Malvern [1969] or Allmendinger et al. [2011]. For large strains it matters whether the strain is described with respect to the undeformed state E -Lagrange-Green strain tensor -or, with respect to the deformed state E * -Almansi-Green strain tensor -: with right Cauchy-Green deformation tensor C and left Cauchy-Green deformation tensor B. The components of finite strain tensors E and E * read [Malvern, 1969]: When the off-diagonal terms of F are small, E and E * reduce to the infinitesimal strain tensor (eq. 15). The finite strain tensors are insensitive to rigid body rotations as C = F T F = QF T QF =F T Q T QF =F TF .

Example: non-proportionality of principal finite strains with length change
The large strain tensors have the advantage that shear-induced extension is accounted for, but they have the drawback that principal strains are not proportional with length change.
For uniaxial extension (a > 1, b = 1, γ > 0, θ = 0), F = Λ, and the Almansi strain tensor becomes: Hence, in the x-direction there is a finite strain (E * xx ) that is not proportional to the stretch: λ = 1 + = a = E * xx + 1. Still, the principal finite strains E * 1 can be related to principal stretch by λ i = 1 + 2E * i [Malvern, 1969], but more intuitive deformation tensors can prove more useful, see the next section.

Large deformation: finite stretch
Instead of measures for strain we will make use of stretch tensors that are a decomposition of the deformation gradient F. Whereas for small deformation, rotation ω and strain can non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International be summed (eq. 14), for large deformation we can no longer assume that shapes stay approximately constant after applying subsequent deformations (i.e. x ≈ X). Thus, instead of a summation as in section 2.5.1, the deformation gradient F can be written as a product of a rotation R (an orthogonal rotation tensor: RR T = I) and a symmetric, stretch tensor with positive principal values [Malvern, 1969]. Thereby R has the form of the rotation matrix Q (equation 13) that applies a rotation with angle θ. Depending on the order of multiplication of rotation and stretching we can write F as: Where U is called the right stretch tensor and V is the left stretch tensor. Figure 4 shows how U describes the shape change in the unrotated state (the original configuration) while V provides the shape change after applying the rotation. Hence, V represents the shape change in the final state, which is why it is our stretch tensor of choice. U and V are related to the Cauchy-Green deformation tensors (eq. 23) by: U and V have the same principal values λ i , and these are related to the principal values µ i of C and B by λ i = √ µ i . The principal axes of U and V are however different, due to the rotation of R, see figure 4. U and V are computed by polar decomposition, for which we apply the Hoger and Carlson [1984] algorithm that provides a closed form solution, see appendix B. Because U and V are decompositions of F, the principal finite stretches are proportional to length changes (as is not the case for finite strain, see previous section). Principal stretches represent the extension in the principal directions, and are in the range [0,∞], with values smaller than 1 denoting shortening and values larger than 1 extension. Because the stretch tensors U and V are symmetric, simple shear will introduce rotations, even in the absence of proper rigid body rotations. R describes the mean rotation, i.e. averaged over the orientations [0, 2π].

Interpretation of the stretch tensor
In 2D, the left-stretch tensor has the form: While the tensor components of V denote extensional (diagonal components) and shearing (off-diagonal components) deformation, the individual tensor components are of limited use for intuitive interpretation of deformation, especially if rotation angle θ (of R) is large. Because V is preceded by rotation R (eq. 28), the deformation described by V does not apply to material that was in the undeformed state oriented along X or Y , but rather to the material 13 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International that is oriented along X or Y after the rotation by R (figure 4). The same arguments hold for U, where the rotation by R follows after the deformation under U.
From a kinematic perspective, we prefer to represent the stretch tensors not in the chosen reference coordinate system, which may be arbitrarily oriented compared to the deformation, but along the principal axes n: where the dyadic product n i ⊗ n i changes the basis, by rotating the reference coordinate system to the directions of the principal axes n. Any second order tensor T can be rotated under a counterclockwise angle to a different base with T = Q T TQ, using rotation matrix Q (eq. 13). Rotating back to the original frame occurs with T = QT Q T . Thus, a more intuitive way of describing the stretch in terms of its eigen values is by changing the basis through a rotation of the coordinate system: where Λ = diag(λ 1 , λ 2 ), and orthogonal rotation matrix Q (equation 13) rotates to the first principal axis n 1 under an angle:

Infinitesimal vs. finite strain
Infinitesimal strain well describes small deformation, which can be understood as deformation where the shape of the deformed state x does not deviate significantly from the reference state X. We have briefly hinted at shortcomings of the infinitesimal strain tensor in section 2.5.1 and will show here in more detail why many common quantities derived from infinitesimal strains are problematic for large deformations. Infinitesimal strains are especially unfit to describe large simple strains or rotations, because the approximation that strain and rotation can be summed no longer holds. A comparison between the the strain tensor and stretch tensor V provides insight when the use of introduces significant errors to the interpretation of strain. Figure 5 makes a comparison in terms of a number of scalar deformation measures: principal strains, maximum shear strain, dilatation and rotation.
As both and V are symmetric tensors, the Mohr circle applies to both, so that we can derive the maximum shear strain as a function of principal strains γ max = 1 − 2 , or principal stretches γ max = λ 1 − λ 2 .
Dilatation, area change, is the product of the principal stretches. In some studies the dilatation is represented by the divergence of the displacement field (i.e. x + y ). This can be non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International polar decomposition of F material axes principal stretches U principal stretches V Figure 4: Sketch of the polar decomposition of the deformation gradient F into an orthogonal rotation R and symmetric stretch tensors V or U. Here we apply the deformation of F to an initially square material element, with side vectors dX. The deformed configuration (right) leads to an element with side vectors dx = FdX. The decomposition can be written as first a stretch U followed by a rotation R, or first a rotation R followed by a stretch V. Principal stretches λ of V provide the stretch in the final, deformed, configuration, whereas stretches λ of U provide the stretch in the original configuration. Material axes provide reference for the reader. We also show strain ellipses in black, that can be constructed using circle coordinates as dX. Because V and U are symmetric tensors, the rotation angle of R is non-zero for simple shear contributions. non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International understood as an approximation to the true dilatation (eq. 4): For small strains only, this reduces to: The last two terms are equal to the first invariant of the strain tensor: The first invariant lacks the 1 2 term in eq. 34, which is why the divergence of the displacement field, or the first strain invariant, is not a good measure for area change in case of large deformation.
As an example we deform material under a simple shear γ, i.e. F ss = Γ (eq. 21). Figure  5 demonstrates that for strains > 0.1, principal strains, rotations and dilatation become increasingly unreliable when using instead of V. For γ > 2 infinitesimal strains even provide physically impossible values for principal strains and dilatation (implying negative lengths or nagative area change). For maximum shear γ max , both infinitesimal and finite tensors provide the same results. Infinitesimal rotation ω increases linearly with shear γ (to the extend of > 2π rotations), while θ correctly converges to π/2, which implies a gradual flattening of material under shear.
non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Figure 5: The limited validity of infinitesimal strain (equation 15) for simple shear with engineering strain γ, i.e. F ss = Γ (eq. 21). Left upper panel: principal strains, 1 , 2 (infinitesimal) and λ 1 ,λ 2 (finite), where 1 does not capture the extension under large simple shears. Also, 2 linearly decreases with increasing γ, past the extensional strain limit of -1. The grey area depicts physically impossible values (section 2.2). Right upper panel: the maximum shear gives the same values for γ max = 1 − 2 (infinitesimal) or γ max = λ 1 − λ 2 (finite). Left lower panel: area change in terms of dilatation (eq. 4). Dilatation using infinitesimal strains dA/dA 0 = ( 1 + 1)( 2 + 1) increasingly deviates for strains > 0.1 from the exact values from dA/dA 0 = λ 1 λ 2 . Moreover the infinitesimal strain can lead to physically impossible values when the dilatation drops below the minimum possible value of dA/dA 0 = 0 (grey area). Lower right panel: mean rotation, using infinitesimal rotation ω (infinitesimal, equation 17) or angle θ from the rotation matrix R (finite, equation 28). Infinitesimal rotation is always proportional to γ for simple shear, while θ converges to π/2, which indicates the flattening of material under simple shear and the alignment of all material orientations with the X-axis. 17 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International

The shape change tensor of choice: V
To quantitatively describe shape changes, the left stretch tensor V, expressed in its principal space (equation 32), is thus most appropriate as it describes length changes as proportional stretches λ. The length changes are provided in the final configuration, and being in the principal space these are independent from the, arbitrary, orientation of the chosen reference frame. The principal values of the finite strain tensors E and E * are not proportional to length change, and are thus a less likely choice for the quantitative description of shape change. Infinitesimal strains are unfit for large strains in the case of rotations or simple shear. In the following we will discuss progressive deformation, after which we will develop qualitative measures of shape change in section 2.7.

Incremental vs. finite strain for progressive deformation
We can write the final deformation F as a product of incremental deformations F i (e.g., Ramberg [1975], Allmendinger et al. [2011]) Because: Even though we can write the final deformation as a product of all previous incremental deformations, the final deformation is still insensitive to the deformation path, as this information is lost in the product of equation 37.

Example: Differences between incremental and cumulative shear strains
Incremental and finite stretches may diverge in terms of the ratio between principal stretches during progressing deformation. This is especially the case for simple shear, but applies to any deformation that contains a shear component. During progressive simple shear (a = 1, b = 1, γ = 0, θ = 0) the incremental deformation gradient takes the form F i = Γ. Considering a constant incremental deformation gradient, the final simple shear deformation gradient F (ss) after n increments becomes: The principal stretches λ, of both U and V, provide us with a description of the stretch that is independent from the coordinate orientation, and can be derived from the principal non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International values µ of C [Malvern, 1969], being solutions of the characteristic equation |C − µI| = (C 11 − µ)(C 22 − µ) − C 12 C 21 = 0. The left Cauchy-Green tensor for simple shear is: The principal stretches λ follow from solving the characteristic equation, using the definition of C (ss) : Here λ 1 is the largest principal stretch. For small deformation λ 1 ≈ −λ 2 , but for increasing n λ 1 will exceed λ 2 in magnitude (see figure 7). Notably, the direction of the largest principal stretch of V, which provides the stretch in the deformed configuration, starts at 45 • for small shears (and thus incremental simple shear) and converges slowly to 0 • for large strains [Allmendinger et al., 2011], showing the gradual flattening of material and stretching in the direction of the X axis.
The dilatation A/A 0 for progressive simple shear (eq. 4): (40) Which confirms that simple shear is area preserving, and that:

Classifying strain type
Different types of incremental tectonic deformation lead to very distinct structures or fault patterns. Therefore it is valuable to infer the type of incremental deformation that has led to the current, finite, deformed state. Especially, if a large fraction of the deformation is localized along discrete faults, we would like to classify the deformation as due to slip on normal, strike-slip, or thrust faults. In our framework we observe all deformation as distributed as PIV averages deformation within a certain spatial window. Our aim is to classify the observed semi-planar deformation as extensional (normal faults), strike-slip (strike-slip faults) and shortening (thrust faults), and the same categories apply to distributed deformation. Here strike-slip is any form of shear, irrespective of rotational components of the deformation.
non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International

Strain type: small deformation
For small, infinitesimal strains we can classify the type of strain by comparing the magnitudes and signs of the largest and smallest principal infinitesimal strain, and the same holds for strain rates. Kreemer et al. [2014] define a strain type measure (˙ 1 +˙ 2 )/max(|˙ 1 |, |˙ 2 |). For strike-slip both principal strain rates will be similar in magnitude but different in sign, leading to a strain type value of 0. For extension and shortening one principal strain rate will be dominant and its value will be positive or negative, leading to strain type values of 1 for extension ( > 1 for extension in two directions). We extend this model into 2 dimensions, as shown by figure 6 where we relate the value and sign of the largest and smallest magnitude principal strains. Using two dimensions instead of a ratio, we aim to develop an intuitive strain measure. We compute the angle φ between the largest and smallest magnitude principal incremental strain: Angle φ thus provides the counterclockwise angle between the two principal values, measured from the max axis. Figure 6 provides a visual interpretation of angle φ, as a function of the two principal strains (equation 42). Angle φ being a discontinuous function, we convert φ into a continuous strain type measure Φ by: Φ thus becomes centered around strike-slip (Φ = 0) and has a range [− π 2 , π 2 ], of which the limits correspond to biaxial shortening vs. biaxial extension.
The strain type Φ is a continuous scale, and within a tectonic context we can expect intermediate, oblique, strain types, notably transtension as an intermediate case between b) and c) having contributions from both strike-slip and extension; and transpression as an intermediate between c) and d) having contributions from both strike-slip and shortening.

Strain type: large deformation
We now consider progressive deformation, leading to large, finite deformation. In considering the evolution of principal strains for progressive, steady deformation, we assume a constant displacement gradient ∂u ∂X . This is most applicable for tectonic deformations, as convergence velocities can be seen as relatively constant over time [DeMets et al., 2010]. This leads to a deformation gradient as function of number of increments n: Figure 7 shows the strain (i.e., λ − 1) evolution for a constant incremental displacement gra-non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International dient, for the cases extension, shortening, strike-slip, and oblique combinations of extension and strike-slip: transtension; and shortening and strike-slip: transpression. Incremental deformation gradients F i are defined using distortion matrices Λ (extension and shortening), Γ (strike-slip) or the product ΛΓ (transtension and transpression). We define simple shear γ along the X-axis, and extensional components along the Y-axis (b = 1; a = 1), see eq. 12.
From the left panel of figure 7, it is clear that for all cases that contain shearing, the largest principal strain λ max − 1 will increasingly become dominant for progressing strain. Thus, the angle φ (dashed lines) is not constant with increasing strain. Therefore we cannot use the same equation as previously, as large shear and extension will have a seemingly similar dominance of the largest principal stretch or strain.
From equation 41 we know that the principal stretches for simple shear/strike-slip relate as λ 2(ss) = 1/λ 1(ss) , which implies that: Therefore, for strike-slip we can also make use of the logarithm of the stretch, and obtain a constant angle φ, irrespective of the magnitude of the finite stretch. This logarithmic stretch is related to Hencky strain, which from the viewpoint of the deformed configuration, is defined as [Xiao et al., 1997]: with principal stretches λ i and the respective eigenvectors n i of V, see also equation 31.
In case of uniaxial extension or shortening, the finite stretch for a constant displacement gradient, or a constant incremental strain, has always one principal stretch equal to one, and thus one zero principal strain. For uniaxial extension or shortening the angle Φ is thus constant, irrespective of the magnitude of deformation. Trivially, the logarithmic principal stretch leads to the same constant Φ. The right panel of figure 7 shows that for a constant incremental displacement gradient and Hencky strains, our strain type measure is constant for strike-slip, as well as extension and shortening. For oblique deformation we obtain a strain type based on the logarithmic Hencky strains that is nearly constant with increasing strain. The change in Φ for oblique deformation is sufficiently small to be able to use the measure for qualitative classification of the strain type. We can thus use the Hencky principal strains to classify the type of strain, in a similar way as for small, incremental strains: Here, i min and i max are the indices of the absolute minimum and absolute maximum principal stretch λ, such that |lnλ(i max )| > |lnλ(i min )|. The continuous strain type Φ is then calculated from φ using eq. 43.

22
non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Left panel: principal finite strains max(|λ − 1|) and min(|λ − 1|) (thinner lines, only displayed if non-zero) and the strain type angle Φ based on finite strain (dashed lines). Here it is clear that Φ is not constant for progressing strain for deformation including shear. Right panel: logarithm of the principal strains lnλ. Here it shows that Φ is constant with increasing deformation for strike-slip, extension and shortening; for transtension and transpression it is nearly constant. Colors correspond to the infinitesimal strains from figure 6.
Another advantage of Hencky strains for large deformation is that lnλ spans a range between [−∞, ∞] instead of [0, ∞]. The Hencky principal strains thus have a more symmetric distribution around 0 (no deformation) compared to λ, where shortening falls in the [0,1] range.
Localized forms of shortening or extension in 3D are caused by shear on faults, but these are shear free deformations in our 2D approach as we lack deformation in the z-direction. Projected on the observed 2D surface pure extension or shortening can then be described (in some rotated reference system) as uniaxial distortions Λ, with b = 1. We expect that in a tectonic context, localized biaxial extensions or shortenings will be rare. For faulting, extension, strike-slip and shortening can be regarded as end-member models [Fossen et al., 1994] and we adopt that approach here. For more distributed deformation we expect that biaxial extension or shortening may occur, for example deformation due to gravitational collapse at topographic gradients. Biaxial deformation that also contains a shear component can not non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International be discriminated from biaxial deformation without shear, when both logarithmic principal stretches have the same sign.

Application: strain type for rotating elements under different types of deformation
As a first test, we apply the strain type classification to individual square elements that we deform using a constant incremental F i , ranging from shortening, via transpression to strikeslip, and subsequently via transtension to extension. We also apply variable rotation to show the independence of the strain type on rigid body rotations. Figure 8 shows the principal stretches and direction of V (left panel) and strain types (right panel) for the final configuration. The first 7 rows show a gradual change from shortening (left) to strike-slip to extension. The 2 last columns contain pure shear and no shape change.
The rows apply varying rotation, with the middle row no rotation, the upper rows clockwise rotation and the lower rows counter-clockwise rotation. Mathematically, the progressive deformation F(n) at epoch n is defined as: The strain types in figure 8 are consistent with the applied deformation gradient, with gradual transitions between shortening through strike-slip to extension. Rotation has, as expected, no effect on the derived strain characterization.

Application: synthetic rotating shear zone
In the following test we define a shear zone (simple shear) that we rotate at the same time, to check whether our methods consistently determine the deformation as strike-slip. We define the shear zone by defining Eulerian displacements δu on fixed, regularly spaced points. First, we define strike-slip displacements d in the local, shear zone co-rotating frame, which are defined at epoch i as: Here ∆u is the final displacement after n epochs, erf the error function, x the coordinate perpendicular to the shear zone, with a range centered around 0, and b is a parameter to scale the width of the shear zone. The displacement vector then becomes: 24 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Figure 8: Synthetic test using individual elements that are deformed using a constant incremental distortion tensorF (equation 48). In the x-direction we vary the distortion parameters, in the y-direction we vary the amount of rotation. Left panel: principal stretches λ from left-stretch tensor V, with the largest principal stretch in black, and the smallest in red. Right panel: strain type Φ based on the logarithmic principal Hencky strains lnλ (equation 46). Distortion parameters used for the first 7 columns, from left to right: a = 1, b = 1 + [−1, −2/3, −1/3, 0, 1/3, 2/3, 1] · 5 · 10 −2 and γ = [0, 1/3, 2/3, 1, 2/3, 1/3, 1] · 5 · 10 −2 . The 8th column with pure shear has a = 1 + 5 · 10 −2 and b = 1/a, γ = 0. The last column has a = 1, b = 1, γ = 0. Rotation angle θ varies linearly between the rows between the values [ 2 3 π, − 2 3 π].
where the first term represents the rotated shear displacements, and the last two terms combined are the rotation contribution. Figure 9 shows the evolution of the principal stretches and principal directions of V, also showing the increasing dominance of λ 1 for progressive deformation. Figure 10 shows the consistently constant strike-slip strain type Φ for the full model and time span, unaffected by rigid body rotations.
non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Figure 9: Synthetic test, representing a counter-clockwise rotating shear zone, with results for 3 epochs: 20, 50 and 100 (final). Principal stretches λ of left-stretch V tensor, with largest λ in black, the smallest in red. As some elements move outside the region for which we describe the Eulerian displacements, we loose part of the model during deformation evolution.

Application: analogue model
In a final step, we analyse incremental displacements derived from successive images of an analogue, laboratory model of a coupled convergent and strike-slip tectonic setting. We have designed a model setup such that different regions of the model experience a wide range of non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International different types of deformation: strike-slip faults (simple shear), normal faulting (extension) and thrust faulting (shortening). In the later stages of the model evolution, a part of the model undergoes a clockwise rigid body rotation, which allows for testing of the strain type analysis to rotations. A horizontal velocity discontinuity approximately half way the model leads to a strike-slip zone, and a fixed indenter is responsible for thrusting in the model domain in front of the indenter, see figure 11. Due to the evolution of topography in front of the indenter, we expect to see normal faulting or distributed extension as well, due to gravitational collapse at topographic gradients. The experimental set-up is relevant for many natural indentation settings, e.g., India or Arabia indentation into Eurasia [Tapponnier et al., 1982, Hubert-Ferrari et al., 2003 or the lateral extrusion of the Alps [Ratschbacher et al., 1991], where shortening in front of an indenter is transferred around the corner into strike-slip deformation along the lateral indenter margin. Our experimental setup merely serves as the proof of concept for our strain analysis; further developed experiments are part of Krstekanić et al. [2020] that compares modeling results to natural cases of curved strike-slip systems at the lateral ends of a rigid indenter.

Analogue model experimental set-up
The laboratory model consists of a single 2 cm thick brittle layer made of dry quartz sand, sieved to grain sizes of 100-300 µm, with density of ρ = 1.500 kg/m3, a coefficient of peak friction P = 0.63 and cohesion of 10-40 Pa [Willingshofer et al., 2018]. Figure 11 shows how we apply a horizontal velocity discontinuity to the quartz sand layer by one mobile and two fixed plastic plates/sheets. The moving plate (green in figure 11) lies below the entire model and is connected to an electric motor that applies a northward pull at a constant velocity of V = 10 cm/h. A narrow cut in the middle of this plastic plate (dashed green line in figure  11) allows for the translation and rotation of the moving plate around a fixed pin (red dot in figure 11). This pin represents the pole of rotation and is fixed to the stable plate (blue in figure  11) at its north-western corner. The third basal plastic sheet, located north of the blue plastic plate, is also fixed (grey in figure 11). Directly north of the blue plate, the grey sheet is placed on top of the moving plate, while in the north-western part of the model it is located below the moving plate (figure 11), which is facilitated by the cut in the moving plate.
Our set-up allows the moving plate to translate into the area north of the pin, imposing a convergent deformation south of the stable region (i.e., blue plate in figure 11) as well as to rotate clockwise. With this configuration, velocity discontinuities are located at two locations: i) at the bottom of the model along the southern and western margin of the stable blue plate and ii) along the contact between moving plate and grey fixed sheet.
Applied displacements in the experiment fall in two distinct phases. In the first phase, the motor pull induces four cm of northward translation, while translation guide bars (left panel of figure 11) prevent rotation. The second phase starts when the moving plate hits the rotation guide bar. At that moment, the translation guide bars are removed, enabling coupled non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International rotational and translational movement of the moving plate (right panel of figure 11). When depressions deeper than 1 mm form, we fill them manually with alternating layers of differently colored sand, representing syn-kinematic sedimentation. We stop the experiment when the total clockwise rotation reaches 4.5 • and the total northward offset along the western margin of the stationary blue plate is 9 cm (configuration as in right panel of figure 11).
The length ratio L * = L model /L nature = 1.5 · 10 −6 , which means that total thickness of 2.0 cm of our model scales to 13.3 km in nature. The strength profile as a function of depth (inset of figure 11) follows the equations by [Brun, 2002], where the strength profile represents the initial conditions.

Tectonic structure analysis and PIV displacement analysis
The PIV analysis is based on top images of the full experiment, taken every 45 seconds. For the analysis we focus on the central part of the experiment (blue outlines of supplemental video S1), which in the images measures 3254 times 1501 pixels. We analyse this area of interest using the software PIVlab [Thielicke and Stamhuis, 2014]. We pre-process images using the PIVlab CLAHE filter and the auto contrast stretch to optimize the contrasts. For the PIV, the first pass uses an interrogation area of 88 px, and a step of 44 px. In subsequent two passes, the interrogation areas are 64, 32 and steps are 32 and 16 px, respectively. We postprocess PIV-results by removing displacements that are larger than 7 times the displacement standard deviation, and manually remove outliers. PIVlab interpolates the resulting missing data.

General model evolution
Initial deformation Figure 12 depicts the principal stretches λ max , λ min and our inferred strain type after 7.5 minutes of model evolution, next to the image at the same time. Principal stretches > 1 imply extensional deformation and are generally shown by the largest principal stretch λ max , while stretches < 1 denote shortening deformation (equation 31). The strain type Φ is based on the relative magnitudes of the logarithmic principal stretches (Hencky strains), see section 2.7. At this point in the model evolution the initial deformation is distributed in broad shear zones: along the N-S velocity discontinuity with transtension bounded at two sides with transpression, rather than by pure strike-slip; indicated in figure 12 by (1). In front of the stable region, along the E-W oriented velocity discontinuity, two zones of shortening evolve. This early structural pattern consists of a top-to-N basal thrust (2) and one top-to-S back-thrust (3).
non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International

Development of faults
At ≈ 1 cm of northward displacement of the moving plate (after 15 minutes), strain localizes into narrow zones of deformation, as figure 13 shows. Along the N-S oriented velocity discontinuity, this early deformation is accommodated by en-echelon NNE-SSW oriented Riedel R-shears (4). Principal stretches λ in figure 13 show NNW-SSE deformation features (5) (P-shear structures), which are characterized by lower strain than the R-shears. Inspection of the photos suggests that the strain levels in the connecting shear zones are too low to result in brittle failure, as no fault structures have formed here. At the same time, the basal thrust (2) and back-thrust (3) become sharper features above the E-W striking boundary of the velocity discontinuity. Around the SW corner of the stable region, this distributed strikeslip deformation (6) connects and transfers into thrusting accommodated by the basal thrust and back-thrust system in which new thrusts start to develop (7). Above the basal thrust, in its immediate hinterland, distributed extensional deformation takes place (8). This extension accommodates the transition between the ramp and the upper flat segments of the basal thrust. With the progress of northward translation, additional en-echelon R-shears form, visible in figure 13 as strain features with the same orientation as (4). The transpressional area at the sides of the shear zone, bounding the Riedel shears, becomes slightly elevated when positive flower structures form.

Merging of strike-slip faults and new back-thrusts
By the end of the northward translation stage (4 cm of moving sheet displacement, at about 30 minutes), all R-shears have become connected to form a continuous strike-slip fault zone where the eastern fault (9) takes up most of the displacement, see figure 14. This process is also visible in the incremental rotation ω in video S2. The previous fault that connected the strike-slip zone to the back-thrusts (6) is well visible in the cumulative deformation, but since the eastern strike-slip fault is now dominant (9), a new fault propagates toward the south (10) to connect with the newest formed back-thrust and thus cuts through the previously formed thrust wedge. The latter fault is short-lived, so its effect on the cumulative deformation is small, but it is well visible in the top-view photo as well as in the incremental strain in video S2. To the SW of the indenter, the back-thrusts become progressively more transpressional (11), as shown by the orange colors in the strain type panel. The distributed extensional deformation on top of the wedge (8) now occurs in large areas above the basal thrust.

Effect of concurrent translation and rotation
After the onset of rotation (starting after 30 minutes) and until the end of the experiment (63.8 minutes), the structural pattern becomes more complex. Figure 15 shows that in the north (i.e., north of the pole of rotation) a transpressional wedge forms, with oblique-slip thrusts (12) bounding the strike-slip deformation zone. In the N-S deformation zone we find the main non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International strike-slip fault (9) that accommodates most of the deformation at the east of the shear zone. Figure 16 clearly shows this combination of strike-slip and thrusting in cross-section a-a' . To the south of the rotation pole, a transtensional basin opens (13) in response to the rotation. There, main dextral-normal faults control the subsidence of the basin, while more distributed oblique-normal slip structures form within the basin ( figure 16).
At the SW corner of the stable region, approximately N-S oriented normal faults (14) with a dextral slip component, accommodate the transition from the highly elevated thrust wedge in the south and the low topography strike-slip zone along the western margin of the stable region. All of these normal faults are connected with the main N-S oriented strikeslip/transtensional fault. At the end of the experiment, the contractional wedge in the south is composed of the basal thrust (2) that accommodates most of the deformation, and numerous faults in the back thrust system (3). The episodic transfer of shortening towards a newly formed back-thrust (15), south of its successor, is visible in video S2 that shows incremental infinitesimal strain and video S4 with strain type. At the same time, distributed extension occurs at the slopes of the wedge (8), see also the visual interpretation from figure  16, cross-section c-c' . We also note three straight and alternating E-W striking features west of the deforming area. These features are not due to real deformation, but rather seem to be a PIV-artifact, caused by interference of the contrasting dark marker lines. intersects the main N-S dextral strike-slip zone at the location of the opened basin, which is bounded by normal faults. Cross-section c-c' intersects the elevated thrusting wedge, showing six back-thrusts, at this location, and the basal thrust that accommodates the largest part of the convergence between the moving domain and the stable region.
non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International

Comparison with visual fault interpretation
To validate that our analysis of the cumulative deformation properly describes the incremental deformation, figure 17 (left panel) overlies the final strain type with the faults inferred from the visual inspection of the model surface and the cross-sections. The location of the backthrusts and basal thrust at both sides of the frontal wedge coincides, as does the the gradual transition from the back-thrusts towards the N-S strike-slip zone. The inferred strike-slip and strike-slip faults correlate well, and the type of obliquity agrees for transpressional and transtensional regions in the N-S shear zone. We also observe distributed extensional deformation on top of the wedge and NW of the corner -with relatively small magnitudes, see the principal stretches in figure 15 -where this is not visible in the images due to absence of localized structures.
As an alternative validation of the strain type, we consider the dilatation field, which in 2D conforms to area change and that provides an additional view on extensional and shortening deformation (whereas it is insensitive to strike-slip). The right panel of figure 17 shows that areas of dilatation > 1 correspond to areas with extensional strain types (or transtensional) and areas with a dilatation < 1 correspond to areas with shortening (or transpressional) strain types.

Time evolution of deformation
We can also analyze the strain type in specific locations as a function of time. This provides a view on when a certain type of deformation was active in a confined material region. Furthermore, it serves as a check how the strain type evolves when different phases of deformation, with possibly different strain types, affect the cumulative strain type. For this purpose the temporal evolution of cumulative strain type and incremental strain type (based on the incremental strain tensor) can be compared. Figures 18 and 18 show the temporal evolution of the principal stretches, dilatation and the strain types for a selection of points in the model.

The basal thrust
Point a (figure 18) is positioned in front of the basal thrust (the actual fault is located where multiple W-E material lines have converged), and the evolution of the logarithm of the principal stretches λ (Hencky strains) shows how after 20 minutes in model time this part of the model is increasingly shortened, ln(λ max ) stays at a small value (i.e. λ max ≈ 1), while ln(λ min ) becomes increasingly negative throughout the rest of the model evolution. From the moment of accumulation of significant strain, the strain type refers to continuous shortening. While the incremental strain type is more noisy, from the moment deformation increases, it indicates shortening in a single direction. The cumulative and incremental strain types are thus in agreement for this location. The dilatation (being < 1) is consistent with shortening as well, non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Figure 17: Final deformation, overlain by faults from visual analysis of fault structures, at the end of the model (64 minutes). Left panel: strain type. Right panel: dilation (relative area change). A dilatation > 1 indicates an area increase (extension/transtension) and a dilatation < 1 indicates an area decrease (shortening/transpression), it is insensitive to pure strike-slip.
as the surface area decreases gradually at this point.

Back thrusts
We focus on two back thrusts, point b focuses on the first back thrust that emerges, point c is located on a later fault. We can discern relatively localized deforming regions, even though non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International our method cannot detect discrete faults. The Hencky strains lnλ show a single direction of shortening, with a clear time delay between the start of the model and the start of shortening. During a later phase (after 20 minutes), the incremental strain type in point b hints at transtensional deformation, as the incremental strain type is generally in between strikeslip and extension (right panel). The dominant strain type in both b and c at the end of the experiment is however shortening, in agreement with the observed thrust faults.

Connection between back thrust and strike-slip zone
From the southernmost back thrust, a thrust fault with oblique slip kinematics (figure 16) connects to the strike-slip zone, and we select here point d. This point is subject to shortening from 10 minutes onward until 20 minutes. From the incremental strain type plot we can see that the deformation becomes more oblique after 20 minutes and gradually becomes strike-slip, which renders the cumulative strain type to oblique shortening (i.e. transpression).

Strike-slip zone
An R-shear, sampled at point e (figure 18), shows strike-slip activity up to 30 minutes, after which the principal stretches stay more or less constant. This specific point shows some signs of transtension, while neighboring points to the west tend to more transpressional deformation. This gradual change in obliquity is well visible in the other direction: points to the east contain a larger fraction of extensional deformation. The strike-slip motion migrates further east after the R-shears connect (see the differences between figures 13 and 14. Point f lies within the main strike-slip fault, and shows most of the time purely strike-slip motion, exemplified by the equal magnitude of the two Hencky strains lnλ. Only at the end the dilatation diverges slightly from 1.

Extension on the thrust wedge
At the edges of the thrust wedge, we find distributed and localized extension. Distributed extension at the northern side of the wedge, where point g initially shows extension as the largest Hencky strain is positive. After 20 minutes the influence of the basal thrust becomes noticeable, and the negative Hencky strain becomes the largest in magnitude. Both dilatation and strain type show this reversal. The incremental strain type, even though it is more noisy, shows this reversal, approximately 10 minutes before the extensional deformation is balanced by shortening deformation in the cumulative strain (right panel). At the transitional region between the high topography wedge and low topography strike-slip region, we find point h, in a region with abundant normal faulting (figure 17). As both Hencky strains are positive, this point is subject to bi-axial extension. The varying orientation of the normal faults in this area, figure 17, also hint at extension in more than one direction. non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Figure 18: Temporal evolution of principal stretches and strain type as well as dilatation of a few selected points: a, along the basal thrust; b, in the first back thrust; c, in a later back thrust; d, along the fault that connects the basal thrusts to the strike-slip shear zone. Upper left panel: strain type (final) and overview of the selected areas. Left row: zoom on the strain type and the selected grid cell (which is outlined in red, neighboring cells outlined in black). Middle row: logarithm of the two principal stretches (Hencky strain) in time. Right row: dilatation (left axis) and strain type (cumulative) and incremental strain type (based on the incremental Lagrange-Green strain, eq. 23) (right axis). 40 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Figure 18: continued. Temporal evolution of principal stretches and strain type as well as dilatation of a few selected points: e, in one of the initial Riedel R-shears; f, along the main strike-slip fault; g, on the wedge close to the basal thrust, where the deformation inverts from extension to shortening (along more or less the same principal direction); h, in the extensional area at the western side of wedge. Left row: zoom on the strain type and the selected grid cell. Middle row: logarithm of the principal stretches (Hencky strain) in time. Right row: dilatation (left axis) and strain type (cumulative) and incremental strain type (based on the incremental Lagrange-Green strain, eq. 23) (right axis). 41 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International 6 Discussion 6.1 Following material during translation, rotations and shape changes The method that we apply to follow the displacement of material through time (e.g. Senatore et al. [2013], Stanier et al. [2016, Boutelier et al. [2019]), building on spatial velocity fields derived from imagery (PIV), makes it possible to study deformation of a small piece of material in an environment that is subject to large translations, rotations or shape changes. Compared to the method of Boutelier et al. [2019] we improve the determination of the displacement gradient by using a bilinear interpretation using shape functions (appendix A) that takes into account the gradual geometry change of the grid, similar as Senatore et al. [2013]. By tracing material points in time we can also determine the cumulative deformation, thereby increasing the signal to noise ratio significantly with respect to the incremental strain. This is shown by our application to PIV-derived velocities of the analogue model of tectonic deformation. The time evolution of the stretch tensor V and the rotation angle θ in supplementary video S3 contains much less noise compared to the incremental strain tensor and the accompanying rotation ω in supplementary video S2. The noise in the incremental deformation seems largely Gaussian in nature, and as a result the cumulative stretch gives much sharper and clearer patterns of deformation. A drawback of using cumulative deformation is that episodical deformation, such as the Riedel shears in the analogue model, or inversions from extension to shortening, can be subdued in the total deformation. Still, it is possible to apply the same method to shorter periods of data to focus on (relatively) homogeneous phases of deformation.

Correct strain tensors for large deformation
We base our deformation analysis on the left-stretch tensor V that gives an exact description of the shape change that is described by the displacement gradient (section 2.5.3). Importantly, V is insensitive to rigid body rotations. The principal values from V (or from rightstretch tensor U) describe the true length changes along the principal axes. Recently, Boutelier et al. [2019]) described PIV-derived finite deformation by U. Contrastingly, the infinitesimal strain tensor (equation 15) that has often been used to describe relatively large deformations (e.g. [Adam et al., 2013, Hoth et al., 2007, Sun et al., 2018, Poppe et al., 2019, Schellart et al., 2019), describes shape changes properly only for small shears and small rotations. Especially, quantities such as principal strains, area change, and rotation from the infinitesimal strain tensor no longer hold for large shear (section 2.5.4). For example, the vorticity ω (eq. 17) applied to large deformations may lead to values that indicate > 2π rotations (figure 5), which is not possible with stationary strike-slip faults. In the absence of additional rigid body rotations, while the rotation value itself implies an incorrect rotation beyond values > π/2, its value can still be proportional to a shear γ. The degree to which infinitesimal strains lead to qualitatively wrong inferences, will thus depend on the setting. Similarly, the use divergence non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International ( x + y ) or dilatation using , to discriminate between areas of area increase or decrease, may lead to misinterpretations for large shear strains (figure 5).

Strain classification
For small strains, the type of strain (e.g. extension, shortening or strike-slip) can be determined from principal strains in a straightforward way [Kreemer et al., 2014]. However, for large strains this distinction between different types of strain can no longer be made based on comparison of principal strains in a similar manner. We introduce the use of Hencky strains (logarithm of finite principal stretch, section 2.7) to characterize strain type. Doing so, we avoid the problem of changing ratios of principal stretches in finite shear deformation (i.e. strike-slip in the horizontal plane), that occurs even for deformation that is constant in time. Hencky strains provide a constant ratio for the two principal strains for extension, strike-slip (simple or pure shear) and shortening, irrespective of the amount of deformation (figure 7). For oblique deformation this ratio only slightly changes with increasing deformation, but for our qualitative measure of strain type we argue this is only a minor issue. Our definition of strain type provides the same qualitative characterization of deformation, as the deformation classification based on visual inspection of fault patterns (figure 17). More so, it provides a fully automated, temporal evolution of deformation magnitude and strain type, and it is able to detect distributed deformation that does not lead to recognizable structures.
When asking the question whether cumulative deformation gives a good sense of episodical deformation related to slip on faults, we see in our analogue experiments that the final strain type agrees with almost all inferred types of faults. Our logarithmic strain measure provides a way of assessing the type of deformation, either distributed, or localized (i.e., faulting), provided that no inversion of deformation takes place. The finite stretch contains no information on the deformation path, so if for example left-lateral strike-slip faulting is followed by rightlateral strike-slip on the same structure, or extension is succeeded by shortening (see point g in figure 18), the cumulative strain type will not be representative for the deformation as occurred during the full observation period. In these instances, it would be preferred to divide the observation period into separate phases (each with its own deformation gradient F) to account for opposite styles of deformation. In a similar way, shortening followed by a phase with strike-slip (see point d in figure 18) is in the finite deformation indistinguishable from concurrent strike-slip and shortening (transpression). Defining different phases of deformation allows for differentiating both options, while still making advantage of the integration of deformation that reduces the noise that may be dominant in incremental deformation.
The interpretation of a strain tensor is non-trivial, as tensor components strongly depend on the orientation of the reference frame (section 2.5.3). Diagonal components and off-diagonal components of the strain tensor are communicative, and vary according to the reference orientation. Often in analogue modeling studies, strain is interpreted by inspecting individual strain tensor components, to derive extensional or shear components (e.g. [Adam et al., non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International 2005, Fischer and Keating, 2005, Hoth et al., 2007, Boutelier and Oncken, 2011, Chen et al., 2015, Le Corvec et al., 2014, Galland et al., 2016, Barcos et al., 2016). This approach is valid only when all deformation zones align with the cartesian reference frame, which limits the use to relatively simple models with straight and perpendicular deformation structures. Contrastingly, principal stretches or strains are unique for a given deformation and have been used to analyze analogue model deformation by Haq and Davis [2009]. Still, principal strains or stretches are inconclusive on shear contributions (in the horizontal plane: strike-slip). Likewise, dilatation can be used as a measure for shortening or extensional deformation [Kettermann et al., 2016], but it is insensitive to shear. In the absence of rigid body rotation, infinitesimal rotation ω (eq. 17) can be used to map small simple shears along strike-slip faults [Le Walter, 2009, Hatem et al., 2017]. However, for large deformations, ω is no longer proportional to the mean rotation of material (figure 5). Even though maximum shear strain γ max will highlight shear zones (e.g. Cruz et al. [2008]), it is just as well sensitive to shortening or extensional deformation. Different from the aforementioned approaches, our method can be applied to any deformation zone orientation, and is independent from the chosen reference frame.

Decomposition of strain
In this study we have developed a qualitative classification of strain, but we do not quantitatively decompose deformation into contributions of shearing and extensional deformation. In recent literature there are multiple efforts to decompose deformation into contributions of shear, extension and rotation. Wang et al. [2019] aim to separate rotations induced by shearing motions and rigid body rotations, based on analysis of the minimum angular velocity in all directions in a small volume, called the Liutex method. Holmedal [2020] applies this model to deformation of solids. The Liutex method however relies on the implicit assumption that shear rotations and rigid body rotations act in the same direction (i.e. the induced vorticity or spin has the same sign). One could imagine a combination of a right-lateral strike-slip fault, inducing a clock-wise rotation within the shear zone, combined with a counter-clockwise rigid body motion. Furthermore, as the Liutex model does not consider extensional motions, even though extensional motions induce internal angular velocities as well [Allmendinger et al., 2011], it is of limited use for decomposing the deformation of solid (compressional) materials.
Other, somewhat older, studies derived decompositions of a 3D deformation gradient into a succesion of a single, volume preserving extension, a simple shear, a rotation and a dilatation [Wang, 1996], or alternatively, a dilation, a rotation and 2 simple shears [Zheng et al., 2000]. As a description of tectonic deformation in the horizontal (2D) plane requires the inclusion of two orthogonal extensions, we cannot use these decompositions to describe tectonic deformation in a physically meaningful way.
A recent, promising decomposition of deformation is the QR decomposition of F that al-44 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International lows for a separation of deformation into simple shear, extension and rotation contributions [Srinivasa, 2012, Freed et al., 2019 (see equations 11 and 12). For 2 dimensions, the extension contribution includes 2 orthogonal components and a single simple shear contribution. In section 2.4 we apply this scheme to describe a general deformation in the direction of the chosen reference cartesian axes. In 2 dimensions, this implies that the decomposed simple shear is always parallel to the x-axis. To describe shear in the horizontal plane on strike-slip faults, we would thus first have to locally align the x-axis of our reference frame with each fault to obtain a meaningful decomposition of deformation. Namely, the QR decomposition is possible in any chosen reference axis orientation [Freed and Srinivasa, 2015]. This implies the development of a spatially and temporally varying definition of the reference axes, and we do not pursue that in this work as it would require the automatic detection of fault structures. Similarly, deriving slip orientations on faults (e.g. Leever et al. [2011]) requires a prior definition of fault orientations. As such this method works only when the spatial derivative of the displacement field is determined perpendicular to a fault.

Outlook
In this study we have made use of an initial grid that is homogeneous throughout the full area of interest. During the subsequent time steps the grid gradually deforms, which leads in general to an inhomogeneous spatial resolution of the derived deformation fields, as extensional or shortening stretches have opposite effects on the grid side lengths. To prevent overly loss of spatial resolution in extensional areas (including shear zones) we have implemented an homogeneous grid refinement with respect to the resolution of the displacement fields as provided by the image correlation techniques. An efficient method to cope with the spatial resolution gradually becoming more inhomogeneous in deforming areas is to use unstructured meshes that are initially inhomogeneous, with increased resolution in deforming regions. This suggests the use of triangular meshes, for which the shape functions can be easily adapted.

Conclusions
We have applied methods to update Lagrangian material displacements using incremental, space-based (Eulerian) displacements from image correlation techniques. We show that the resulting displacement fields can be used to compute high-resolution, time-dependent, 2D shape changes -described as strains or stretches -even for large deformations.
In a similar way as for small strains, it is possible to use principal stretches to classify the type of strain that has led to finite deformation. For this purpose, we introduce a novel use of logarithmic principal stretches, also known as Hencky strains, to qualitatively describe 2D plane deformation. Hencky strains have the desired property that the relative lengths between the two principal strains do not change during constant incremental deformation. non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International Tectonic deformation, i.e. shortening, strike-slip or extension, and their transpressional and transtensional transitions, can be adequately described using our developed qualitative strain type measure. Due to the temporal and spatial continuous description, this measure can be used to inspect the temporal evolution of shape changes, even in the absence of discrete structures.

Acknowledgements
We have made use of color maps from ColorBrewer2.0, by Cynthia A. Brewer and Scientific Colourmaps [Crameri, 2018]. The analogue model experiment has been conducted in the framework of the PhD project of Nemanja Krstekanić. 53 non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International For the validity of the bilinear shape functions the requirement holds that all angles between the element sides need to be less than 180 • , which in practice means elements should be sufficiently small to prevent occurrence of angles larger than 180 • . Our elements are isoparametric, which means that the shape functions not only interpolate the coordinates, but serve also as interpolation functions for other nodal quantities, such as displacement.

CReDiT
We have the displacement values u n available at each of the four nodal points, and we can interpolate these displacements anywhere within the quadrilateral element using our shape functions N: u(s, t) = u n N(s, t) where u = u v T and nodal displacements are: Our main purpose here is to determine ∂u ∂X , that we can write as: We can thus compute the displacement gradient as the product of the displacement gradient with respect to the local coordinates, and the inverse of the gradient of the global coordinates X to the local coordinates ∂X ∂s , which is often called the jacobian J. Using equation 53 we compute the elements of ∂u ∂s at the centroid of the element (s = 0, t = 0): non-peer reviewed EarthArXiv preprint submitted for review to Geophysical Journal International

B Polar decomposition of the deformation gradient in 2D
Analytic solutions exist for the decomposition of the deformation gradient F into the rotation tensor R and the right stretch tensor U as well as the left stretch tensor V [Hoger and Carlson, 1984]. As the polar decomposition holds: Knowledge on U and U −1 leads to V and R (recall that RR T = I) by: Hoger and Carlson [1984] provide expressions for U in closed form as well as its inverse, without the need to determine the tensor square root of the right Cauchy-Green tensor C, as U = √ C.
Firstly the principal values µ of C: µ 1 = C 11 + C 22 + (C 11 − C 22 ) 2 + 4C 2 12 2 (68) From these the invariants of C and U can be calculated: By inserting 68 into 70 we can further simplify the invariants and bypass the computation of principal values: Then the right stretch tensor is: and its inverse: Which we can insert in equation 67 again to obtain V and R. The rotation angle θ represented by R can be calculated using the four-quadrant inverse tangent: