Enhanced Digital Image Correlation Analysis of Ruptures with Enforced Traction Continuity Conditions Across Interfaces

: Accurate measurements of displacements around opening or interfacial shear cracks (shear ruptures) are challenging when digital image correlation (DIC) is used to quantify strain and stress fields around such cracks. This study presents an algorithm to locally adjust the displacements computed by DIC near frictional interfaces of shear ruptures, in order for the local stress fields to satisfy the continuity of tractions across the interface. In the algorithm, the stresses near the interface are extrapolated by local polynomials that are constructed using a constrained inversion. This inversion is such that the traction continuity (TC) conditions are satisfied at the interface while simultaneously matching the displacements produced by the DIC solution at the pixels closest to the center of the subset, where the DIC fields are more accurate. We apply the algorithm to displacement fields of experimental shear ruptures obtained using a local DIC approach and show that the algorithm produces the desired continuous traction field across the interface. The experimental data are also used to examine the sensitivity of the algorithm against different geometrical parameters related to construction of the polynomials in order to avoid artifacts in the stress field.


Introduction
Understanding the dynamics of shear rupture of interfaces is important for fields ranging from failure of composites and bimaterial structures, to earthquakes, car brakes, and even to pistons of internal combustion engines. Recent advances in high-speed camera technologies have enabled the development of an experimental setup that combines ultra-high-speed photography with digital image correlation (DIC) to quantify the behavior of propagating dynamic ruptures in the laboratory [1][2][3], including full-field measurements of displacements, velocities, strains, and stresses ( Figure 1).
DIC is an optical technique used to determine the full-field deformation between a reference and a deformed image [4][5][6]. The correlation algorithms use either global or local approaches to compute the desired field quantities [7,8]. In global approaches, the pattern matching is performed over the entire domain of analysis, typically using the finite element method [9,10], which allows enforcing compatibility of the displacement field. In local approaches, the image matching is performed with small windows, or "subsets", separated by a distance, referred to as "step", which can be less than half a subset size, (i.e., subsets can overlap) [4,11,12]. A number of recent investigations have focused on evaluating the accuracy and resolution of DIC [13,14]. The full-field measurements enabled by DIC have been used in a wide range of applications [4,6]. However, standard DIC techniques cannot capture displacement discontinuities associated with cracks or ruptures, as they assume a continuous displacement field [4]. Various approaches have been proposed to analyze discontinuous displacement fields [15][16][17][18][19][20][21][22], but they generally entail the application of constraints with a theoretical interpretation of the fields. This limits the implementation of such approaches to cases where the theoretical assumptions are valid.
In a previous study on dynamic ruptures conducted using an experimental setup similar to that employed here, it was possible to capture the static displacements associated with dynamic ruptures, after these ruptures traversed a part of a fault, because a low-frame-rate camera with high resolution (i.e., 4 megapixels) was employed [23]. In that study, the discontinuous displacement field was accurately mapped by simply having subsets over the interface as the high resolution and low noise of the images allowed for the subset size to be relatively small compared to the image size, so that smearing of the displacement fields caused by having subsets across the interface was minimal. Generally, local DIC approaches that provide the solution up to half a subset from the interface may also be sufficient when tracking the quasi-static propagation and opening of tensile crack faces and the associated strain field, using a large field of views and high-resolution cameras. However, these correlation approaches are not possible when the subset size is large compared to the field of view (i.e., when using low-resolution cameras), such as our ultra-high-speed camera (which has a resolution of 250 × 400 pixels 2 ) or high-resolution cameras, but with comparatively large subset sizes, situations typically arising in the study of dynamic problems [3]. To resolve the displacement discontinuities on the interface associated with the propagation of dynamic shear ruptures, the studies in [1][2][3] used the commercial DIC software Vic-2D (Correlated Solutions, Inc.) to separately correlate the domains above and below the interface. While standard local DIC methods calculate displacements up to half a subset away from the interface, the "fill boundary" algorithm of Vic-2D employed in [1][2][3] uses affine transformation functions to extrapolate the displacements from the center of the subset up to the interface. This enables the use of large subsets in each side of the interface, which are essential to overcoming the noise associated with ultra-high-speed photography and the analyzed rapid ruptures. An important limitation of this approach, however, is that the extrapolated displacements are less accurate than those obtained directly by the actual correlation, especially at pixels near the interface, which have the largest distance from the subset center. This leads to errors in the strains and stresses near the interface, which are obtained from the spatial derivatives of the displacement fields. In addition, the non-coupled correlation of the domains above and below the interface results in non-physical discontinuities in the tractions across the interface.
This study proposes a simple and fast method to supplement the DIC solution with a postprocessing algorithm that enforces the continuity of normal and shear tractions across the interface. The algorithm uses a constrained inversion to construct local 2-D polynomials that satisfy traction continuity conditions at the interface while matching the displacements of pixels closer to the center of the subset, where the DIC solution is most accurate. The polynomials are then used to calculate all stress components near the interface. In Section 2, we describe the experimental setup, including the laboratory experiment, ultra-high-speed diagnostics, the digital image correlation approach, and the post-processing procedure to turn the displacement fields into strain and stress fields. In Section 3, we present the method to enforce traction continuity along the interface and study the effects of the parameters involved in this method on the full-field stresses, displacements and particle velocities. Implications for the analysis of friction using the stress fields produced by this approach and conclusions are given in Sections 4 and 5, respectively. Figure 1. Schematics of the laboratory rupture experiment. Dynamic shear ruptures evolved spontaneously along the frictional interface, inclined at an angle α, of two Homalite plates under a compressional load P. Ruptures were initiated by the small burst of a NiCr wire placed across the interface and connected to a capacitor bank. The white light produced by a flash source was reflected by the specimen's surface and captured by a low-noise high-speed camera, typically at 1-2 million frames/s. The portion of the specimen to be imaged, the field of view, was coated by a flat white paint and decorated by a characteristic speckle pattern. Next, the textured images were processed by digital image correlation (DIC) algorithms to produce a temporal sequence of full-field displacement maps. The displacement fields were then post-processed to produce velocity, strain, strain rate, and stress maps.

Monitoring Dynamic Shear Ruptures in the Laboratory
The algorithm developed in this paper is a post-processing algorithm that is designed to enforce the continuity of stresses measured near frictional interfaces during experimental dynamic ruptures. This section briefly summarizes the experimental setup, diagnostics, image analysis, and numerical methods that are used to obtain the displacement, strain, and stress fields before using the algorithm. This setup is the evolution of the Caltech " Laboratory Earthquake Setup " developed by Rosakis and his co-workers over a span of 15 years [1][2][3][24][25][26][27][28][29]. A more detailed description of the current form of the setup is given in [3].
The presentation focuses on the analysis of experimental dynamic ruptures, but the algorithm can be used for a wider range of frictional experiments that involve image analysis. The same methodology can also be used for the analysis of opening, shear, and mixed-mode cracks at both coherent and incoherent (frictional) interfaces separating both similar and dissimilar solids. Such interfaces are common in nature (e.g., faults), and are also important to a variety of engineering problems involving composites and bi-materials.

The Laboratory Setup
The laboratory setup is designed to study the dynamics of shear ruptures propagating along preexisting inclined frictional interfaces via full-field measurements of displacements, velocities, strains, and stresses associated with the rupture (Figure 1). Two Homalite-100 quadrilateral plates with a frictional interface inclined at an angle α are loaded under uniaxial compression P (Figure 1), NiCr wire resulting in initial shear and normal stresses on the interface of τ0 = Psin(α)cos(α) and σ0 = Pcos 2 (α), respectively. The rupture is nucleated by a local pressure release provided by a rapid expansion of a NiCr wire filament due to an electrical discharge of Cordin 640 high-voltage capacitor (Cordin, Salt Lake City, UT, USA). A key aspect of the setup is that the low shear modulus of Homalite enables production of well-developed dynamic ruptures in samples of tens of centimeters.
Once the rupture initiates, a target area coated with a random black-speckle pattern is monitored using Shimadzu HPV-X ultra-high-speed camera system (Shimadzu, Kyoto, Japan), capable of recording up to 10 million frames per second, and Cordin 605 high-speed white-light source system with two light heads (Cordin, Salt Lake City, UT, USA) ( Figure 1). In the experiment reported here, the camera recorded a sequence of 128 images of the patterns distorted by the propagating rupture, with a resolution of 400 × 250 pixels 2 , at a temporal sampling of 2 million frames/second and an exposure time of 200 ns. Moreover, the HPV-X camera was equipped with a range of prime telephoto lenses, which allowed it to monitor fields of view of different sizes [3].

Digital Image Correlation to obtain Displacement Fields
We employed the local digital image correlation (DIC) software Vic-2D (Correlation Solutions Inc.) to analyze the sequence of images acquired with the HPV-X high-speed camera, and to produce evolving displacement maps, computed with respect to a selected reference configuration. In local approaches, the correlation is performed on local "subsets" separated by a distance, referred to as "step", which can be less than half the subset size (i.e., subsets can overlap) [4,11,12]. The 2D-DIC algorithms provide the two in-plane displacement components at each subset center. A parametric study performed in [3] showed that a subset size of 41 pixels and step size of 1 are needed to accurately resolve the spatial and temporal features of dynamic ruptures. In order to capture the discontinuous displacement field across the interface, the correlation is performed separately for the domains above and below the interface. While standard local DIC approaches are able to produce the displacement map up to half a subset away from the interface, the "fill boundary" algorithm of Vic-2D uses affine transformation functions to extrapolate the displacements from the center of the subset up to the interface.

Post-Processing of the Displacement Fields
Prior to the computation of strain fields, we filtered high-frequency noise from the displacement fields using a non-local-means (NL-means) filter [23,30,31]. This filter enables the displacement fields to be smoothed, while maintaining the original signal pattern. An example of full-field displacement maps of a laboratory rupture obtained using high-speed photography, DIC analysis, and a non-local filter is given in Figure 1. To facilitate comparisons, we report the same experimental rupture discussed in [1,3], produced with an applied vertical load P = 23 MPa and inclination angle α = 29°, and monitored with a small field of view. This was a super-shear rupture propagating at a speed of 2200 m/s. Super-shear ruptures are interfacial ruptures whose tip propagates at speeds greater than the shear wave speed of the surrounding solid, as first observed in [24,26,32,33]. Note that the fullfield maps are cropped and their size is slightly smaller than the reported field of view size. These displacement fields are used throughout the study to test the capability of new algorithms to enforce traction continuity (TC) at the interface.
Strains are computed from the filtered displacement fields using finite difference approximation. We employ a frame system x1-x2 parallel to the interface and use the central difference scheme for pixels away from the boundaries [23]: , , 1, , 1, , 2 , , 1 2 where , , and , , are the interface-parallel and interface-normal displacement components, respectively, for pixel , and frame k, and s is the step size. Note that the displacements in Equations (1-3) are expressed in pixels. Second-order backward and forward difference schemes are used to compute strains at the pixels immediately above and below the interface, respectively [1,3].
The stress changes with respect to the reference configuration (before rupture) are computed from the strain fields using the standard plane-stress linear elastic constitutive equations: where E is the Young's modulus and ν is Poisson's ratio. The actual stresses are given by: , , , where , , , , and , are the pre-stresses at the reference configuration. Since Homalite-100 is a strain-rate sensitive material at the strain rate levels developed during these dynamic ruptures, we use the dynamic values of the elastic constants to compute the stress changes [1,3]. An example of full-field maps of normal and shear stresses calculated from the displacement fields shown in Figure 1 are given in Figures 2a and b, respectively. The figures clearly show that the experimental technique enables capturing the spatial and temporal features of the ruptures. In particular, the more compressional and dilatational regions surrounding the crack are clearly imaged, as expected for shear ruptures (see [1,3] for more examples and physical insight). However, a significant limitation is that in the "fill boundary" algorithm of Vic-2D the domains above and below the interface are correlated and extrapolated towards the interface independently, without applying continuity constraints along the interface. Therefore, tractions are discontinuous at the interface, especially the normal component, where at some locations there is a difference of about 9 MPa between the normal tractions above and below the interface. Note that the increase of normal tractions towards the interface is the physical consequence of the shear rupture; however, close to the interface the tractions need to evolve towards being continuous. In the following section we present a post-processing algorithm to resolve this issue.

A Post-Processing Algorithm to Enforce Traction Continuity along the Interface
We present a method to supplement the DIC solution with a post-processing algorithm in which the displacements are modified in the proximity of the interface, such that the normal and shear tractions are continuous across the interface. The main idea is (i) to keep the displacement fields obtained with the local DIC up (or close) to the center of the subset touching the interface, and (ii) to use a portion of those fields to construct a polynomial extrapolation for the displacement so that stresses computed from the extrapolated displacements satisfy traction continuity conditions at the interface. The formulation presented here is for a coordinate frame x1-x2, such that x1 is parallel to the interface. However, the method can easily be extended to any other frame. The procedure is described here for a given column of pixels at a given frame, and was implemented for all the columns and frames.

Traction Continuity Conditions
The condition for continuity of normal traction across the interface is given by  where the superscripts ( 0+ ) and ( 0-) denote the values immediately above and below the interface, respectively, and are the initial stresses. Assuming a linear elastic response for the materials above and below the interface, the condition in Equation (6) can be written in terms of strains as where the superscripts ( + ) and ( -) denote the material properties above and below the interface, respectively. For small strains the condition can be expressed in terms of displacement gradients as Similarly, the continuity of shear traction is expressed via and , where is the shear modulus.

Approximating the Displacements with Local Polynomials
In order to modify the displacements near the interface for the j th column of pixels, we approximate the displacement fields , , , and with local polynomials in terms of both spatial coordinates. The polynomials are defined over rectangles of pixels above and below the interface, which are centered at the j th column of pixels ( Figure 3). The pixels within each rectangle are divided into a group of , pixels that are closer to the interface, and a group of , pixels that are closer to the center of the subset. To represent the spatial variations in stresses, we approximate the displacement fields with cubic polynomials with n = 10 coefficients as , Note that and should be chosen small enough to avoid smoothing local patterns in the displacement fields.

Inverting for the Polynomial Coefficients
To obtain the polynomial coefficients, we aim to minimize the misfit between the observed displacements u , u , u , and u , produced by the DIC analysis, and the displacements predicted by the polynomials u , u , u , and u at the pixels closer to the center of the subset. Recall that the DIC solution is computed up to half a subset away from the interface and extrapolated to the interface by the "fill Boundary" algorithm of Vic-2D. Since we expect such extrapolation to be reliable close to the center of the subset, but not close to the interface, we can use some of the extrapolated values by taking where are the sensitivity sub-matrices above and below the interface, respectively, are the assemblies of polynomial coefficients, and , is the number of pixels within the group closer to the center of the subset.
The polynomials must also fulfill the constraints for continuity of tractions across the interface. The condition for continuity of normal tractions is enforced by and can be assembled into the following matrix form: where, for example, Similarly, the condition for continuity of shear stresses is enforced by Different methods can be used to solve the system of equations in Equation (12) together with the constraints in Equations (14) and (15) for the polynomial coefficients; here we use Matlab (The MathWorks, Natick, MA, USA) function lsqlin. Once the coefficients are obtained, the displacements at the pixels of the central column ( pixels above and pixels below the interface) are replaced by those predicted by the corresponding polynomials, and the procedure is repeated for all columns and images. Note that the system of equations to be solved is small because the local polynomials should represent the local variations in displacements and stresses well. Thus, the computation time was small, and the modification of all pixel columns within an image was done within 2-3 s for the images analyzed here (400 pixels along the interface). Figure 4 shows the interface-parallel (a) and normal (b) displacements at the center of the field of view (j = 190) for the rupture of Figure 1. Both the displacements obtained by Vic-2D (with the "fill boundary" algorithm) and displacements modified by TC enforcement are shown. The maximum difference in displacements between the two solutions was smaller than 0.001 pixels, which is equivalent to 0.5% for u1 and 1.5% for u2, but it was enough to enforce the continuity of traction. When the gradients in displacements in the direction parallel to the interface are large, as in the case of the rupture analyzed in this study, direct calculation of strains and stresses using the modified displacement fields may not result in completely continuous tractions at the interface. The calculation of strains and stresses for pixels in column j involves the displacements in columns j−1 and j + 1, as it requires displacement gradients with respect to x1 (Equations 1-4). When enforcing TC, these gradients are computed using the polynomial centered at column j. Once the displacement fields have been updated using the procedure described above, strains can be computed from the updated displacement fields, using finite difference schemes, as discussed in Section 2.3. However, the finite differences would be using the displacements of columns j−1 and j + 1, which are obtained from other polynomials than that used to enforce TC, those centered in columns j−1 and j + 1. As such, the spatial derivatives of the modified displacement field with respect to x1 would be slightly different than those obtained during the enforcement of TC, using the polynomial centered at column j. Although the differences between polynomials centered at two neighboring columns are small, they may be enough to cause some discontinuities in tractions at the interface. To ensure the continuity of tractions, we computed strains and stresses within the 2 band, where displacements are updated directly as the algorithm goes over each column of pixels (i.e., we calculated the displacement derivatives for column j with respect to x1 using the displacements in columns j−1 and j + 1 obtained from the polynomial centered at column j. Full-field maps of stresses modified by the algorithm using nx1 = 7, nx2,inter = 12, and nx2,cent = 11 (Figure 2 c,d) demonstrated the capability of this approach to generate more realistic and continuous stress fields near the interface. Interestingly, because of the symmetry of the problem, the enforcement of TC led to normal stresses that are close to the background value (∼17.6 MPa) at the interface.

The Effects of the Geometrical Parameters of the Polynomial
It is important to choose appropriate values for the geometrical parameters nx1, nx2,inter, and nx2,cent in order not to introduce unwanted field distortions. For instance, a choice of a small value of nx2,inter together with a large value of nx2,cent or vice versa can lead to stress fields that are continuous on the interface, but have sharp discontinuities (Figure 2 e,f) or spurious oscillations (Figure 2 g,h) in the bulk. In this section, we study the effects of different parameter combinations in more details, aiming to find the parameter sets that avoid these distortions in the stress field. We ran the algorithm 324 times, varying the geometrical parameters in the following ranges: 3 ≤ nx1 ≤ 9, 4 ≤ nx2,inter ≤ 20, and 5 ≤ nx2,cent ≤ 21. Note that for nx1 ≥ 13 the system was over constrained. Overall, we found that parameters in ranges of 8 ≤ nx2,inter ≤ 14 (1/5 to 1/3 of the subset size) and 7 ≤ nx2,cent ≤ 15 gave similar results, independent of the value of nx1.
We found that the stress fields obtained with the algorithm were insensitive to the value of nx1. For example, there was no observable difference between full-field maps of stresses obtained with nx1 = 3 and nx1 = 9 ( Figure 5 a-d). To further examine the effect of nx1, we plot the normal and shear stresses vs position along two lines perpendicular to the interface at x1 = 170 and 220 pixels (Figure 5e-h). For both lines there was a "jump" larger than 5.5 MPa in the normal stress at the interface when the calculation was done with the "fill boundary" algorithm of Vic-2D without TC enforcement. The TC algorithm enforced the normal stress to be continuous across the interface, with a normal stress at the interface that was about the average of the normal stresses at the pixels immediately above and below the interface when TC was not enforced. The results were not sensitive to nx1. A small effect of nx1 on the shear stress was observed for the line located at x1 = 220 pixels (Figure 5h), with a difference of 0.12 MPa between the results for nx1 = 3 and 9. The shear stress at the interface when TC was enforced were smaller than the shear stresses at the pixels immediately above and below the interface when TC was not enforced. The stresses calculated with the "fill boundary" algorithm of Vic-2D without the traction continuity algorithm are also shown, where the black diamonds correspond to the centers of the subsets above and below the interface. In all the cases shown in this figure, nx2,inter = 10 and nx2,cent = 11.
As shown in Figure 2 e-h, a large difference between the parameters nx2,inter and nx2,cent resulted in artifacts in the stress field. On the one hand, taking nx2,inter = 4 and nx2,cent = 21 (with nx1 = 7) resulted in discontinuities at the boundary between the Vic-2D solution and polynomial extrapolation (  x 2 (pixels) x 1 (pixels) n x1 = 9 n x2,inter = 10 n x2,cent = 11 x 1 (pixels) n x1 = 9 n x2,inter = 10 n x2,cent = 11 4.7 9.3 14 e,f), nx2,inter + nx2,cent pixels away from the interface due to a poor fit between the stresses obtained by the polynomial and the original stresses at this region. On the other hand, taking nx2,inter = 20 and nx2,cent = 5 (with nx1 = 7) resulted in spurious oscillations (Figure 2 g,h). Some effects of nx2,inter were observed also when nx2,cent was fixed at a moderate value (e.g., nx2,cent = 11), especially for the shear stress. A small value of nx2,inter (e.g., nx2,inter = 6) led to large gradients in the normal stresses near the interface, as well as to small discontinuities at the pixels located nx2 pixels above and below the interface (Figure 6 a,b). A large value of nx2,inter (e.g., nx2,inter = 16) resulted in an artificial increase or decrease of stresses at localized regions around the interface compared with the domains above and below (Figures 6 c,d).
For example, the normal stress decreased locally at the interface at x1 = 370 pixels (Figure 6c), while the shear stress increased locally at x1 = 170 and 370 pixels (Figure 6d). Some issues are illustrated further in plots of the normal and shear stresses vs position along two lines perpendicular to the interface at x1 = 170 and 220 pixels (Figures 6 e-h). The value of normal stress on the interface was similar for nx2,inter ≤ 14, but diverged for higher values, with a difference of about 1 MPa for nx2,inter = 20. For nx2,inter ≤ 6, the TC enforcement resulted in a poor agreement between the stresses obtained by the polynomial and the original stresses also at the pixels closer to the center of the subset. The shear stress at the interface also diverged, with an increase of about 1 MPa for nx2,inter = 16, at x1 = 170 pixels ( Figure 6f). The effect of nx2,cent was generally smaller than that of nx2,inter. Full-field maps of stresses obtained with nx2,cent = 7 and nx2,inter = 10 (Figures 7 a,b), were similar to those in Figures 2 c,d. However, large value of nx2,cent = 19 resulted in small discontinuities in normal stress at the pixels located nx2 pixels above and below the interface (Figures 7 c,d). Both at x1 = 170 and 220 pixels there was a very small effect of nx2,cent on the normal stress at the interface (Figures 7 e,g). However, for nx2,cent ≥ 17, the fit between the stresses obtained by the polynomial and the original stresses was poor also for the pixels closer to the center of the subset. Similarly to nx2,inter, the shear stress at the interface may significantly have differed for large values of nx2,cent (Figure 7f).

Implications for Friction Analysis
The experimental setup developed in [1,3], which combines ultra-high-speed photography with digital image correlation, enables us to observe displacements, velocities, and stresses close to the interface. For the field of view shown in Figure 1, these quantities can be measured up to 25 μm (0.5 pixels) from the interface using the "fill boundary" algorithm of Vic-2D. This enables studying the evolution of friction, defined by the ratio / , at any point along the interface, as well as its dependence on variables such as slip (relative displacement across the interface) and slip rate [1]. As a first approach to enforce traction continuity on the interface, the tractions and were calculated in [1] by averaging the stresses at the pixels immediately above and below the interface. Because of the anti-symmetric and symmetric patterns in the stress changes ∆ and ∆ , respectively, the normal stress on the interface was inferred to be nearly constant, as expected during the propagation of a dynamic shear rupture along a planar interface, and the shear stress had a smoother evolution than the individual components and above and below the interface. The averaging enabled us to study friction, despite large differences between the normal stresses immediately above and below the interface [1], but it is important to examine whether those average values actually correspond to the tractions at the interface computed with the continuity condition.
In general, the averaged normal stresses (without TC enforcement) agreed with the normal stresses computed with TC enforcement for the pixels immediately above and below the interface (Figure 8a), with some small deviations. Note that, because we enforced the continuity at the interface, the stresses obtained with TC enforcement for the pixels located 0.5 pixels above and below the interface did not completely overlap, with a maximum difference of 0.3 MPa. Because of its symmetry with respect to the interface, the discontinuity in shear stress across the interface for the original Vic-2d results (Figure 8b) was significantly smaller than in the normal stress. Overall, there was a good agreement between the shear tractions obtained by the averaging procedure and those obtained with TC enforcement, again, with some minor differences. Interestingly, the differences were consistently smaller in the cohesive zone (the zone where the shear traction drops from a maximum to a residual value), with a maximum difference of about 1 MPa at x1 = 9 mm. This is also shown locally in Figure 5h. Because the stresses were computed from displacement gradients, small changes in the displacements were enough to enforce the traction continuity on the interface, and there were only minor differences between the interface-parallel displacements and velocities obtained with TC enforcement and those obtained without TC enforcement (Figures 8 c,d). The small differences in stresses resulted in some differences in the evolution of local friction with slip and slip rate (Figures 8 e,f). However, the important characteristics of the local friction, such as the peak, residual levels, and weakening distance, were similar. Figure 8. (a, b) normal and shear stresses along the interface. The black curves represent the stresses ±0.5 pixels from the interface obtained without traction continuity (TC) enforcement, while the green curve represents the average of the two. The blue curves represent the stresses obtained with TC enforcement ±0.5 pixels from the interface. (c, d) Interface-parallel displacements and velocities ±0.5 pixels from the interface, obtained with (blue) and without (black) TC enforcement. (e, f) Resolved friction vs slip and slip rate at x1 = 8.5 mm, obtained with (blue) and without (black) TC enforcement. Both curves represent the average of the values immediately above and below interface. The polynomial geometrical parameters are nx1 = 7, nx2,inter = 10, and nx2,cent = 11.

Conclusions
In this work, the local DIC solution is supplemented with a fast post-processing algorithm to enforce the continuity of tractions across the interfaces of shear ruptures. This procedure allows us to obtain more physically meaningful stress fields near the interface, which is important when DIC is applied to study the dynamics of laboratory frictional ruptures. In the algorithm, the stresses near the interface are calculated from local polynomials that are constructed using a constrained inversion. This inversion is such that the resulting displacements match the displacements of pixels closer to the center of the subset, where the DIC solution is most accurate, while the resulting stresses satisfy traction continuity conditions at the interface.
We applied the algorithm to the displacement fields of a laboratory shear rupture obtained with a local DIC procedure, to show that the algorithm indeed produces the desired continuous traction fields across the interface. A sensitivity study provided a constraint on the parameters nx2,inter, nx2,cent, and nx1 involved in the construction of the polynomials, such that undesired artifacts in the stress fields were eliminated. We found that parameter ranges of 8 ≤ nx2,inter ≤ 14 (1/5 to 1/3 of the subset size) and 7 ≤ nx2,cent ≤ 15 gave similar results, independent of the value of nx1. Relatively minor changes in displacement fields around the interface can produce non-negligible gradient changes, resulting in some non-uniqueness regarding the exact stress evolution towards the interface, even if the traction continuity is enforced. Future progress can be made by combining the results inferred by the analysis presented in this work with dynamic rupture modeling that matches the observed full fields. The dynamic solutions can then be interrogated for the spatial dependencies of stress fields next to the interface.
The averaging procedure for stresses above and below the interface, employed in the previous study [1] to study friction, works well to mimic the traction continuity across the interface due to the special nature of this rupture problem that exhibits certain symmetries across the interface. However, the averaging procedure may not work as well in cases where the symmetries are disrupted, as in the case of ruptures approaching a free surface [27]. We will examine such cases in future studies.