## 3 3D Kinematics of the Boundary Displacement Field  
To characterize the 3D kinematic properties of displacement fields, it is a prerequisite to define the 3D basic unit, which kinematic operators can apply. Figure 4 gives a schematic illustration of a 3D basic unit that is defined for the present kinematic analysis. In this study, a series of auxiliary origins along the axial direction of the specimen are assumed, for which the radial and tangential displacements are assumed to be zero; the axial displacements are evaluated as the averaged vertical displacements sampled along the boundary at each specimen height. This assumption is inspired and supported by the experimental observations of [PERSON] et al. [7]: according to the X-ray tomography results that the deformed sand specimen consists of a rigid cone located along the axis of the specimen, and several shear planes extending outwards from the cone to the specimen boundary, which implies that the radial and tangential displacements along the axis are negligible compared to the displacements along the vertical direction. The defined auxiliary origins are then linked with the boundary coordinates used in the DIC analysis, to form the basic unit for kinematic analysis in cylindrical coordinates. In a preceding paper by [PERSON] and [PERSON] [3], the authors have provided the first-order kinematic operators in cylindrical coordinates. Specifically, the first-order kinematics comprised of _gradient_, _divergence_, and _curl_ fields that are produced by applying the _nabla_ operator to the vector field \(U\)[22]. The _gradient_ is a tensor product between the  
Figure 4: Schematic illustration of the 3D basic unit for kinematic analysis in cylindrical coordinates. Figure 3: Incremental displacement fields of test 092903b between an axial strain of 3.0% and 9.0% (post-peak of the stress–strain curve). The first row from left to right are displacements along horizontal (\(u\)), vertical (\(v\)), and out-of-plane (\(w\)) directions. The second row describes the same displacement field, but decomposed into radial (\(r\)), tangential (\(t\)), and axial (\(v\)) directions in cylindrical coordinates. _nablaola_ operator and the 3D vector field \(U\) (Equation (2)), which includes nine derivative components representing translational or rotational deformations of a local area. The sum of three diagonal terms is the _divergence_ of the vector field \(U\) (Equation (3)), which indicates the volumetric sink (negative quantities) or source (positive quantities) of a local area. The _curl_ contains three terms \(\left(curl\,U\right)_{\rho}\), \(\left(curl\,U\right)_{\phi}\), and \(\left(curl\,U\right)_{y}\) representing local rotational deformations along the \(\hat{\rho}\), \(\hat{\phi}\), and \(\hat{y}\) axes, respectively (Equation (4)):  
\[grad\,U=\
abla\otimes U= \begin{bmatrix}F_{11}&F_{12}&F_{13}\\ F_{21}&F_{22}&F_{23}\\ F_{31}&F_{32}&F_{33}\end{bmatrix} = \begin{bmatrix}\frac{\partial U_{\rho}}{ \partial\rho}&\frac{1}{\rho}\left(\frac{\partial U_{\rho}}{\partial\rho}-U_{ \phi}\right)&\frac{\partial U_{\rho}}{\partial\psi}\\ \frac{\partial U_{\phi}}{\partial\rho}&\frac{1}{\rho}\left(\frac{\partial U _{\phi}}{\partial\psi}+U_{\rho}\right)&\frac{\partial U_{\phi}}{\partial\psi }\\ \frac{\partial U_{\phi}}{\partial\rho}&\frac{1}{\rho}\frac{\partial U_{\phi}}{ \partial\psi}&\frac{\partial U_{\phi}}{\partial\psi}\end{bmatrix}