## 2 Experimental Method
### 3D Digital Image Correlation  
The digital image correlation (DIC) method is a non-destructive, optical method for measuring displacement fields over the surface of a deforming body. The object tracked by DIC analysis is a cluster of image pixels forming a subset. When studying material evolved from initial state to deformed state, DIC measures the affine behaviors (i.e., translation, rotation, and/or straining) of overlapping subsets between the reference and target images. Due to the innate heterogeneity of sand colors, the distinct greyscale patterns of each subset can be directly recognized through a matching algorithm without any manual intervention. In this work, two digital cameras were obliquely set up in front of the specimen to construct a 3D scene, which is similar to the way human eyes acquire the object position [1]. The calibration of lenses involves capture-synchronous images of a standard grid oriented at different angles. From the resulting stereo images, key camera parameters (e.g., position, orientations, focal length, lenses distortions, etc.) are calibrated. After confirming the object's 3D spatial information, 3D displacements can be computed by comparing overlapping pixel subsets between reference and target images as described above. The software VIC-3D, by Correlated Solutions, Inc. (Correlated Solutions: Irmo, SC, USA), was used to extract 3D displacement fields from the stereo images. The search for best-match subsets was implemented following normalized cross-correlation algorithm [32]. Once matching subsets and corresponding displacement vectors were identified, a continuous displacement field was linearly interpolated to evaluate subset translation, rotational, and straining properties on a continuum domain. Details of the piece-wise integration of cumulative displacement fields can be found in [25]. The stereo images were  
Figure 1: Global material responses: (**a**) triaxial stress-strain curves; (**b**) axial strain-volumetric strain curves. taken every 15 s, corresponding to 0.05% of axial strain through the course of compression. An increment of 0.2% of axial strain, however, was considered sufficient to represent the evolution of localization effects [25; 27]; thus, incremental displacements were updated every 4 images. The measurement accuracy was assessed by comparing averaged vertical displacement through DIC data with global readings by linear variable differential transformer (LVDT) transducers, and the difference was found approximately along the vertical direction. Similar accuracy has been observed in the horizontal direction [1; 33]. To illustrate the components of a typical displacement field over the specimen's surface, Figure 2 presents displacement fields of test 092903b between axial strain 3.2% and 7.0% using the Lagrangian description resulting from 3D DIC analysis. The first row shows displacement vectors decomposed into horizontal (\(u\)), vertical (\(v\)), and out-of-plane (\(w\)) directions, respectively, based on a Euclidean coordinate system. The horizontal displacement, as anticipated, related to the development of an expansion band, which was concentrated in the middle of the specimen and made the sand move towards opposite directions (left negative, right positive) if viewed on the horizontal plane. The vertical displacement field (middle sub-plot in the same row) showed approximately zero displacement at the top (since loading was provided from the bottom) suggesting that the specimen was compressed approximately as three separate \"moving blocks\" that possessed different deforming rates [29]. The bulging effect of the specimen is depicted in the plot of the out-of-plane displacement field (third sub-plot), in which the maximum value was shown at the center of specimen. The second row in Figure 2 is the same displacement vector, but decomposed under cylindrical coordinates, yielding radial (\(r\)), tangential (\(t\)), and vertical (\(v\)) displacement fields, respectively (left to right).