Calibration of the Discrete Element Method and Modeling of Shortening Experiments

The discrete element method (DEM) is becoming widely accepted as an effective method for addressing tectonic problems in granular materials. It is capable of reproducing structures observed in the analog model (AM). However, the previous experiments also pointed to variability among DEM models and AMs in the number of fault zones, their dip angle and spacing, and the evolution of the surface slope of a thrust wedge. The accuracy of the DEM depends on the input parameter values, so the calibration of the discrete element method is very important. Microscopic properties of particles and macroscopic properties of loose quartz sand were calibrated through a series of repose angle and biaxial tests. Furthermore, an AM was constructed to simulate the evolution of the thrust wedge to compare with DEM results. DEM and AM results indicate an encouraging overall agreement in model evolution. Based on a new automated wedge quantification method, DEM results were systematically compared with AM results on the number of fault zones, their dip angle and spacing, the evolution of the surface slope of a thrust wedge, and other parameters. Our study provides a necessary comparison between commonly applied modeling approaches, which is important for more confidently applying these methods to understand real fold and thrust belt systems.


INTRODUCTION
Fold and thrust belts are a series of mountainous foothills adjacent to an orogenic belt, which forms due to contractional tectonics. As the total shortening increases in a fold and thrust belt, the belt propagates into its foreland. The frontal parts of these fold and thrust belt systems have commonly been interpreted to be critically tapered Coulomb wedges (Chapple, 1978), which were quantified by Davis et al. (1983), Dahlen (1990), and others. The modeling experiments, such as an analog model (Suppe, 2007;Wu et al., 2010;Wu and McClay, 2011;Schreurs et al., 2016;Sun et al., 2016) and discrete element method (Buiter, 2012;Morgan, 2015;Buiter et al., 2016;Li, 2019), are widely used for quantitative comparison with natural examples of fold and thrust belt systems.
The discrete element method (DEM) is a non-continuum method and is capable of reproducing structures observed in the analog model (AM) and has been applied to the study of geological problems (Saltzer and Pollard, 1992;Hardy and Finch, 2006;Hardy, 2008;Hardy et al., 2009;Yin et al., 2009;Buiter, 2012;Dean et al., 2013;Morgan, 2015;Botter et al., 2016;Morgan and Bangs, 2017;Smart and Ferrill, 2018;Meng and Hodgetts, 2019). Buiter et al. (2016) presented the quantitative results from three brittle thrust wedge experiments of eleven numerical codes, which use finite element, finite difference, boundary element, and distinct element techniques. The continuum methods, e.g., finite element model, finite difference model, and boundary element model, can simulate discontinuities to an extent either replacing the discontinuities with the material of a different rheology or through special treatments of the discontinuity nodes. However, they are not suitable for studying emergent behavior and probe the evolution of fracture systems, which is a consequence of microscopic processes (Gray et al., 2014;Mora et al., 2015). Particularly, the results of Buiter et al. (2016) showed the SDEM, named the stress-based discrete element method, is capable of reproducing structures observed in the analog sandbox experiments without the need for the ad hoc calibration routines normally associated with the conventional DEM. In contrast to the conventional DEM, the friction properties of an SDEM particle system are in agreement with the Mohr-Coulomb constitutive model with friction angles specified on a particle level (Egholm et al., 2007). However, the experiments also pointed to variability among models in the number of shear zones, their dip angle and spacing, and the evolution of the surface slope of thrust wedges. And they lacked quantitative comparison between the SDEM and corresponding AM results. Hardy et al. (2009) presented that a series of numerical simulations at the km scale based on the DEM replicated the similar deformation seen in the AM at the cm scale. The accuracy of the DEM largely depends on the input parameter values, so the calibration of the discrete element method is very important (Coetzee, 2016;Coetzee, 2017). Generally, the particle properties and the emergent bulk material properties are typically assessed through the repose angle tests and biaxial tests (Botter et al., 2014;Morgan, 2015), so that particle aggregates can produce a realistic fault geometry and finite strain field.
In Discrete Element Method, the basic methodology of DEM is outlined. In Calibrating the Microscopic Parameters of Particles, the microscopic properties of particles and macroscopic properties of loose quartz sand are calibrated through a series of repose angle tests and biaxial tests. The experimental setup, the model construction technique, and the material properties are described in Shortening Experiment. Furthermore, an AM using loose quartz sand was constructed to simulate the evolution of the thrust wedge in order to compare with the results of the DEM at the same scales. In Comparisons With Analog Model and Discrete Element Method, an accurate method for measuring the surface slope, width, and height of the thrust wedge is proposed based on the mesh. The results of discrete element simulations were compared with scaled analog (sandbox) models through the determination of their qualitative (visual) and quantitative (e.g., measurements of the surface slope and dip angle of faults) similarities and differences, which is beneficial to the Earth Sciences community.

DISCRETE ELEMENT METHOD
The DEM has been applied to the study of geological problems in recent decades. In standard DEM simulations, the perfect discoid or spherical shape of the particles grossly exaggerates rolling when compared to, e.g., nonspherical and rough quartz sand grains. To compensate for this effect, it is important to incorporate rolling resistivity. This approach also leads to a higher efficient computation. Furthermore, the geometries that occur in the DEM will be characterized by forward thrust propagation and by back thrusting incorporating rolling resistivity. It is similar to what we observe in a typical analog model experiment. The DEM described here is a variant of the lattice solid model, which was developed to provide a basis to study the physics of rocks and the nonlinear dynamics of earthquakes at the beginning (Mora and Place, 1993;Mora and Place, 1994;Mora and Place, 1998;Mora and Place, 1999;Place et al., 2002). Generally, the lattice solid model does not include particle resistance to rolling, i.e., the particle spins are initially set to zero and fixed in the whole simulation, which can enhance interlocking force between particles. The particle shapes of the loose quartz sand are various, and the occlusion between them is strong. Although the discoid particles, i.e., twodimensional disks, are used in our simulation, the particle spins are fixed in order to enhance interlocking force between particles. Furthermore, the fault-fold structures shown in the simulation based on the lattice solid model are similar to those observed in the AM (Hardy et al., 2009). The geologic body is simplified into an assemblage of discoid elements, which obey Newton's equations of motion and can move under the action of forces that are generated by interacting with pairs of elastic springs. A full, detailed description of the theory behind this modeling approach and its application to geological problems are given by Place et al. (2002) and Hardy et al. (2009).
In the simulations presented here, inter-particle bonding is not used, i.e., all springs are not tensile. A full detailed description of the bond is given by Hardy et al. (2009). The behavior of the elements assumes that the particles interact through a repulsive force in which the magnitude of the normal force, F n , is given by where K n is the normal stiffness of the contact and U n is the normal relative displacement. In addition to treating the normal force between particles, we also calculate the shear force, F s , as a result of displacement perpendicular to the vector connecting the particles' centroids, U s . The magnitude of the shear force is limited to be less than or equal to the maximum shear force, μF n , allowed by Coulomb friction, where K s is the shear stiffness of the contact, F n is the normal force at a contact, and μ is the inter-particle friction coefficient. If a contact is "lost" between two touching elements (i.e., they separate), then all the forces between the elements are set to zero. The value of normal stiffness of the contact is set to be equal to the shear stiffness of the contact, and they were not distinguished in the following discussion, i.e., k K n K s . An artificial viscosity (F v ) is added to damp the reflected waves from the boundary of the particle and to avoid buildup of kinetic energy in the closed system (Place et al., 2002;Liu C. et al., 2013). The viscous force is proportional to the particle velocities and is given by F v ηv p , where η is the artificial viscosity and v p is the particle velocity. The time step Δt is a constant value. It is chosen to ensure the stability and accuracy of the numerical simulation, particularly the integration. It is determined on the basis of the maximum stiffness of the contact, k, and the particle with the smallest particle mass, m. The relationship often used is of the form Δt f · 2m/k √ , where f is the safety factor of the time step (Itasca Consulting, 2008).

CALIBRATING THE MICROSCOPIC PARAMETERS OF PARTICLES
The accuracy of the DEM depends largely on the microscopic properties of particles. The microscopic properties of particles and the macroscopic properties of loose quartz sand were calibrated through a series of repose angle tests and biaxial tests, in order to reproduce the specified mechanical properties of loose quartz. To do this, the static repose angle of the material was measured by the repose angle tests, where initially the assemblage is inside a rigid box, and then one side of the box is removed, as shown in Figure 1. The microscopic parameters of particles used in the DEM are given in Table 1. When the particles' friction coefficient, μ, changes from 0.0 to 0.5, the static repose angle of the sample is proportional to the particles' friction coefficient ( Figure 2).
From repose angle tests, Hardy et al. (2009) had determined that their inter-particle coefficients of friction, 0.2, 0.1, and 0.05, correspond to bulk coefficients of friction of ∼0.55, 0.36, and 0.18, respectively (angles of ∼30°, 20°, and 10°), values that allow for comparison with analog modeling studies and that can be compared to natural examples (McClay and Whitehouse, 2004). Additional ledge simulation was used to measure the angle of repose for our model (Oger et al., 1998;Botter et al., 2014). Finally, the friction coefficients of the particle in our model were 0.30 corresponding to the repose angle of ca. 40°, which is also consistent with that of the quartz sand used in our laboratory ( Figure 1).
The dimensions of biaxial samples are 1 × 2 cm ( Figure 3D). The Mohr-Coulomb failure envelope of the material is derived from biaxial tests under various confining pressures (50, 100, 200, and 300 Pa) ( Figure 4). The initial biaxial sample is created by radius expansion (Itasca Consulting, 2008). First, as shown in Figure 3A, 1, 600 particles are created inside a rigid box with smaller diameters (0.1, 0.2, 0.25, and 0.3 mm). Then, their diameters are multiplied by 2  Table 1.
Note: The particle packing consists of four particle sizes, with diameters and quantity ratio of 0.2 mm, 0.4 mm, 0.5 mm, and 0.6 mm and 2:2:1:1, respectively. d, largest particle diameter. ρ, particle density. g, gravitational acceleration. f, safety factor of the time step. k, stiffness of the contact. μ, friction coefficient. η, dynamic viscosity. υ, velocity of the mobile wall.
FIGURE 2 | The static repose angle, θ, of the sample is proportional to the particles' friction coefficient, μ.
Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 636512 ( Figure 3B). Finally, the assembly is cycled to equilibrium. The particles outside the dashed line are deleted ( Figure 3C) to create new sidewalls, and the initial biaxial sample is obtained ( Figure 3D). When the particles' friction coefficient, μ, is 0.0, 0.1, 0.2, 0.3, 0.4, and 0.5, the samples' friction angle, φ, and cohesion values, Co, are shown in Figure 5. They both increase with the increase in the particles' friction coefficient, μ. Klinkmüller et al. (2016) reported the material properties of 26 granular analog materials used in 14 analog modeling laboratories. The peak friction angle and cohesion values measured from ring-shear tests are shown in Figure 6, and the five-pointed star is measured from biaxial tests by the DEM ( Figure 4B). The bulk cohesion value of the DEM is ca. 67.4 Pa. by and large, consistent with the cohesion values of the quartz sand used in analog experiments ( Figure 6).
The DEM for geological applications typically uses larger particle sizes and fewer particles (Saltzer and Pollard, 1992;Burbidge and Braun, 2002;Strayer and Suppe, 2002;Finch et al., 2003;Yamada, 2004;Naylor et al., 2005;Benesh et al., 2007;Hardy et al., 2009;Miyakawa et al., 2010;Hardy, 2011;Vidal-Royo et al., 2011;Hardy, 2015;Morgan, 2015), so that it can model structures at a larger spatial scale applicable to common geological problems. The particle sizes in our models were selected to directly compare against the quartz sand. The quartz sand, with a grain size of 0.2-0.4 mm and the spontaneous stacking density of about 1,500 kg/m 3 , was used in the AM experiment, and its deformation obeys the Mohr-Coulomb failure criterion (Lohrmann et al., 2003;Panien et al., 2006). The average bulk material density of particle assemblage is ca. 1,500 kg m −3 . The assembly with four different particle diameters (0.2, 0.4, 0.5, and 0.6 mm) was used in the DEM. Furthermore, the sizes of the elements used in our DEM (0.2-0.6 mm) are very close to the real size of the quartz sand (0.2-0.4 mm) in our AM. Although the size of disks could be modeled more accurately by using smaller particles, the representation chosen was thought to be accurate enough and yet computationally efficient. The particle aggregate can basically represent the physical characteristics of analog materials used in analog modeling laboratories. A full discussion of parameter testing can be found in the following references: Place et al., 2002;Hardy and Finch, 2006;Hardy et al., 2009;Hardy, 2011;Vidal-Royo et al., 2011;Hardy, 2013;Botter et al., 2014. SHORTENING EXPERIMENT

DEM Setup
Approximately 120,000 particles with four different particle diameters (0.2, 0.4, 0.5, and 0.6 mm, their quantity ratio is 2:2: 1:1) were initialized by randomly generating particles within a 60.0 cm long and 10.0 cm high domain. Then, they were allowed to reach a static equilibrium state under the action of gravity. The particles whose height exceeded 2.5 cm would be deleted, and they were permitted to relax to a stable equilibrium state by settling under gravity again, i.e., all particles in the assembly have come to rest and their positions change only insignificantly, and the kinetic energy and the gravitational potential energy are approximately constant. The resulting particle assembly is 60.0 cm long and 2.5 cm high, and the number of particles is 101, 599. Microscopic parameters in Table 1 were chosen as the particle properties for the simulations. Deformation is illustrated in relation to 10 initially flat-lying layers that are pure for visualization purposes, which allows a much higher resolution of structures to be achieved and observed compared to previous studies of thrust wedges (Burbidge and Braun, 2002;Naylor et al., 2005;Hardy et al., 2009;Morgan, 2015). The computing restrictions limit the number of particles. The particle numbers (ca. 120,000) are carefully chosen to balance the resolution of the calculation domain and the computing restrictions. The left (mobile), right, and basal walls of the model have the same coefficient of friction as particles. Our simulation was run for 129,048 time steps with output of the assembly every 8,065 time steps equivalent to 1 cm (Δt × v × 8,065) of shortening along the base wall, which took ca. 3 h to complete by VBOX Li et al., 2018;Li, 2019), which was compiled with GCC (Richard M. Stallman and the GCC Developer Community, 2012) at O2 levels of optimization and The geometric evolution of the DEM simulation is shown in Figure 7.

AM Setup
The AM experiments, which have a long history in tectonic modeling, form an excellent basis for testing how well the DEM reproduces structures that are relatively well understood. A typical shortening experiment based on the AM was carried out using an apparatus with two glass sidewalls, i.e., a fixed backstop and a mobile wall, all of which were installed on a horizontal Plexiglas base plate. Shortening was induced by the mobile wall, which was connected to a motor-driven piston. The internal dimension of the apparatus was 100 × 30 × 30 cm (length, height, width), while the size of our initial model was 60 × 2.5 × 30 cm. The modeling materials were sieved into the box from a height of 30 cm. Ten initially flat-lying sand layers were sieved on the apparatus (Figure 7). Before sieving, the inside of the apparatus was carefully cleaned using an alcoholic solution and dried thoroughly to minimize sidewall friction. The shortening velocity of the mobile wall was set as 0.004 mm/s. During the experiment, sequential deformation was monitored by digital cameras. Photographs were taken from the side at intervals of 2 min. The structural evolution of the AM is shown in Figure 7 with plots of the geometry after 2 cm, 4 cm, 8 cm, 12 cm, and 16 cm of shortening, and fault zones can readily be identified on the geometry plots. Accretion of the material was in sequence and  can be described in two stages: an initial stage of rapid thickening involving closely spaced thrusts ( Figures 7A-C) followed by a stage of cyclical growth with alternations between wedge lengthening and wedge thickening involving widely spaced thrusts ( Figures 7D,E). General evolution of fold and thrust belts has been the focus of previous studies (Mulugeta, 1988;Liu et al., 1992;Burbidge and Braun, 2002;Schreurs et al., 2006;Bonnet et al., 2007;Bose et al., 2009;Bigi et al., 2010;Wu and McClay, 2011;Schreurs et al., 2016;Sun et al., 2016); the comparison of our model results with previous studies demonstrates that the fold and thrust belts were characterized by forward thrust propagation and by back thrusting. In particular, a shortening experiment performed by Schreurs et al. (2016) had good agreement with our experiment, as shown in Figure 8.

COMPARISONS WITH ANALOG MODEL AND DISCRETE ELEMENT METHOD
Deformation Fields and Strain Analysis Schreurs et al. (2006) and Schreurs et al. (2016) presented the results of an analog comparison study with many participating modeling laboratories, which shows that when one modeler repeats the same experiment, quantitative model results still show variability. In the AM, the deformation should be quantified using image analysis [e.g., particle image velocimetry (Adam et al., 2005;Leever et al., 2011;Boutelier and Cruden, 2017)], laser scanning (miniature laser altimetry) (Liu Z. et al., 2013), and volumetric scanning using computerized X-ray tomography (Suppe, 1983;Boyer et al., 1989;Suppe and Medwedeff, 1990;Colletta et al., 1991;Uehara and Takahashi, 2014;Zulauf et al., 2016). Laboratories that do not have access to tomographic techniques can only monitor the outside boundary of models, with the observation of internal structures limited to cutting the model at the end of the experiment.  The DEM is easy and exactly reproducible under the self-same initial and boundary conditions, allows great flexibility in the choice of geometries and material properties, and can quantify a large range of parameters directly, such as stress, strain, and velocity. The DEM gives detailed information on the displacement trajectories of all the particles in the system. Using this information, a representative, or average, strain for the overall system of particles or a subdomain of the system can be determined. The cumulative displacements of the particles were derived for each new particle configuration and resolved into horizontal and vertical components as shown by Morgan (2015). A triangulation-based nearest neighbor interpolation algorithm (MathWorks, 2015) was used to grid ca. 0.58 mm spacing, which is triple the average radius of all the particles, and a continuous surface was fit to each component. To analyze the simulation processes, involving strains are not small compared to unity, and the theory of finite strain must be used (Oertel, 1996;Means, 1976). The distortional strain, i.e., strain-induced distortion, was used to quantify the results for the DEM and was calculated according to Morgan (2015), and it can be quantified as the second invariant of the deviatoric finite strain tensor.
The structural evolutions of the models are shown in Figure 7 with plots of deformation and cumulative distortional strain after 2 cm, 4 cm, 8 cm, 12 cm, and 16 cm of shortening. We can discriminate more elaborate expressions of fault damage from deformation fields by strain analysis. As shortening went on, the deformation zones in the AM and DEM expanded bidirectionally, i.e., toward both the foreland and the hinterland near the backstop. The hinterland of the AM near the backstop was over a relatively flat terrain; on the other hand, there was not an obvious flat hinterland in the DEM. The inner parts of the proto wedges were progressively involved in deformation by sequential back thrusting identified from distortional strain fields in Figure 7. The deformation zones expanded bidirectionally, toward both the foreland and the hinterland near the backstop, and the "X" shear zones (blue and red zones in the map of distortional strain) come into being ( Figure 7A). Synchronously, the deformation fronts migrated forward, and new thrusts developed in the foreland. Backward thrusts (blue zone of distortional strain in Figure 7) that occurred in the DEM showed relatively little displacement, and the slip was predominantly along forward thrusts.

Fault Interpretation
A back thrust fault and nine forward thrusts in the AM ( Figure 9A) and a main back thrust fault and seven forward thrusts in the DEM ( Figure 9B) were interpreted, involving closely spaced thrusts near the mobile wall (between F1-F6 and F7 for the AM in Figure 9A, between F1-F4 and F5 for the DEM in Figure 9B) and widely spaced thrusts near the foreland (between F7, F8, and F9 for the AM in Figure 9A, between F5, F6, and F7 for the DEM in Figure 9B). Forward thrusts accommodate most of the shortening in both the AM and the DEM, which was consistent with previous discrete element simulations (Morgan, 2015). F0, F3, F5, F7, F8, F9 of the AM and F0, F3, F4, F5, F6, F7 of the DEM were essentially in the same position, as shown in Figure 10A, respectively. Meanwhile, the spaces between F8 and F9 for the AM and F6 and F7 for the DEM were nearly the same ( Figure 10A). Dip angles of the forward thrusts increased from the foreland (ca. 20°) to the mobile wall (ca. 85°), as shown in Figure 10B. This phenomenon is very common among Piedmont tectonic belts, such as the Longmen Shan fold-thrust belt in Central China (Jia et al., 2006;Hubbard and Shaw, 2009) and the South Tianshan Mountain range in the Kuqa region of Northwest China (Wen et al., 2017). Meanwhile, the fault geometry changed from a "line" (e.g., F6-F9 in Figure 9A) to "S" (e.g., F1-F4 in Figure 9A). From the distortional strain field of the DEM ( Figure 9C), linear faults near the foreland and "S" faults near the mobile wall could be identified easily as the strain was very concentrated to help us identify the shape of the faults. Both the steepening of faults Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 636512 toward the hinterland and the change in shape from linear to more "S" shaped reflect the process of break-forward imbrication and fault-related folding of structures to the hinterland in a break-forward sequence of thrusts (Shaw et al., 2005).

Wedge Sensitivities
In order to allow for a comparison of AM and DEM results in a more quantitative manner, the following properties: surface slope, wedge width, and wedge height, were measured. The surface slope, wedge width, and wedge height of the AM were obtained from the cross section shown in Figure 7. So far, there seems to be no general agreement on how to measure thrust wedge surface slopes (Buiter, 2012). The initial surface slope angles are difficult to determine because there is only one thrust . For two or more thrusts, the surface slopes have been determined by drawing the enveloping surface as  Figure 9A and DEM in Figure 9B. The position of a fault was defined as the intersection of this fault and the stratigraphic line between the above two stratums, i.e., the position of F9 in Figure 9A is 29.2 cm and the position of F7 in Figure 9B is 28.0 cm. In particular, the position of the back thrusting, i.e., F0 in Figure 9A and F0 in Figure 9B, was defined as the intersection of this fault and the stratigraphic line between the below two stratums. (B) Dip angle of faults obtained from the AM in Figure 9A and DEM in Figure 9B. As the shape of some faults, e.g., F1 in Figure 9A, is "S." The dip angle of a fault was defined as the dip angle of the segment across the intersections of this fault and the stratigraphic line between the above two stratums and the below two stratums. Therefore, the dip angle of F1 in Figure 9A is ca. 78°, i.e., the dip angle of the yellow line. Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 636512 8 shown in Figure 11A. Oscillations in surface slope angles occur just before or after the formation of a new thrust, depending on the degree to which a new thrust is incorporated into the wedge. Figure 11A shows how they have been determined in our paper. The experimental evolution of the AM and DEM was quantitatively analyzed by those parameters in the following sections. Inevitably, the measurements of surface slopes may be influenced by the measurer and subject to small differences in interpretation . Therefore, in a general way, the values of surface slopes of the AM were completely measured by three people, in the same manner, and subsequently averaged.
Here, an automated wedge quantification method is proposed based on the mesh and illustrated by a given example.
The contractional wedges consist of the top, toe, surface, surface slope, width, and height, as shown in Figure 11A. To find the top point, the coordinates of the point in the surface should be found. In the discrete element model, the object of the study is discrete particles. Here, the points of the surface are searched based on the mesh, and the toe of the slope is found by the farthest distance method. The surface slope, wedge width, and wedge height have been determined by the following steps.

1) Surface search method based on mesh.
First, parallel to the Y-axis, the mesh is divided ( Figure 12A). The mesh width is generally two to three times the maximum particle radius. The highest particle of each mesh was found. Their IDs were stored in an array from left to right. The surfaces were formed by connecting the particles in the array in turn, i.e., the yellow and blue particles in Figure 12A.
2) The farthest distance method searches the toe of the slope.
Among the particles forming the surface, the position of the particle with the maximum y value, T m , is the slope top. The particle with the maximum x value, T n , is the termination point of the surface extension. Among the blue particles between T m and T n , the point T j with the maximum distance from the line segment T m T n is defined as the slope toe ( Figure 12A).

3) Surface slope calculation method.
After the slope toe T j is obtained, the particles between the slope top T m and the slope toe T j can be determined (purple and green particles in Figure 12B). Yellow particles in the middle 80% region of the line segment T m T j are taken to fit a straight line, and the angle between the straight line and the horizontal plane is defined as the surface slope, as shown in Figure 12B.
Then, the width and height of the wedge can be easily calculated from the top and toe of the wedge. The sequential deformation was monitored by digital cameras in the AM. A high-quality binary image of the AM can be obtained based on the digital image processing methods, including binarization and noise removal . The high-quality binary image consists of a rectangular array of pixels. Those pixels can be regarded as discrete square particles. Then, the above method can still be used to acquire the actual surface of the AM.
The same as the surface slopes of the DEM from 4 to 6 cm, the surface slopes of the AM have an obvious fluctuation from 3 to 5 cm, as shown in Figure 11B, based on the new wedge quantification method we presented. The surface slopes of the DEM change slightly after 6 cm, and the surface slopes of the AM change slightly after 5 cm. In the first stage, the thickness of the wedge increased slowly, and the surface slope gradually increased, i.e., when model shortening was from 3 to 4 cm in the AM and from 4 to 5 cm in the DEM ( Figure 11B). In the second stage, i.e., when the model was pushed from 4 to 5 cm in the AM and from 5 to 6 cm in the DEM, a new fault began to grow and wedge lengthening played a dominant role in this stage, which made the surface slopes get reduced ( Figure 11B).
With respect to the surface slope, the wedge width and height of the AM and DEM are more consistent, as shown in Figures  11C,D. The fluctuation of the maximum wedge width in the whole experiments of the AM and DEM indicates that the experiments entered a stage of cyclical growth with alternations between wedge lengthening and wedge thickening involving widely spaced thrusts ( Figure 11C). In the maximum height curves, the maximum wedge height of the AM and DEM with a highly consistent degree increased linearly with shortening, basically ( Figure 11D). The growth of the wedge height is relatively stable, in the whole experiments of the AM and DEM. The results indicate an encouraging overall agreement in model evolution.

CONCLUSION
An automated wedge quantification method is proposed based on the mesh. The results show the width, height, and surface slope of the thrust wedge can be precisely quantified using the proposed Frontiers in Earth Science | www.frontiersin.org May 2021 | Volume 9 | Article 636512 method. It is helpful for accurately measuring the width, height, and surface slope of the thrust wedge for studying the overall dynamic evolution of the numerical and analog models. Although the manner in which the AM and DEM accommodate shortening leads to a similar style of deformation, it is also clear that variations exist among them. Similarities and differences between AM and DEM results are listed as follows.
Similarities: Shortening is accommodated by in-sequence forward propagation of thrusts (Figure 7). The shape of thrusts is consistent, "S" for the smaller thrusts near the mobile wall and "line" for the bigger thrusts near the foreland (Figure 7). The distance between a newly formed thrust and the previously formed thrust is highly variable for small thrusts near the mobile wall (F1∼F6 of the AM and F1∼F5 of the DEM in Figure 9). The positions and the dip angle of forward thrusts near the foreland (F7, F8, F9 of the AM and F5, F6, F7 of the DEM) show variations in Figure 10. The wedge height and width of the AM and DEM increase linearly with shortening, and they are fit very well.
Differences: The number of thrusts that have been formed at a specified amount of displacement is variable, e.g., the number of thrusts in the AM near the mobile wall is more than that in the DEM (Figure 9). The surface slopes remain in the stable field for critical taper theory. The wedge steepens until the critical taper is attained, i.e., the angle between the base and the surface of the wedge reaches a critical value. When the critical taper is attained, the wedge is considered to be stable. The surface slopes of the AM achieve a critical taper soon, but they moderately increase to a critical taper in the DEM ( Figure 11B).
The microscopic properties of particles and macroscopic properties of loose quartz sand were calibrated through a series of repose angle and biaxial tests. Then, a twodimensional DEM with high resolution and an AM experiment were constructed to simulate the evolution of the thrust wedge. Overall dynamic evolution of the DEM and AM is highly consistent, in spite of the difficulty of achieving an exact representation of the analog conditions with the DEM. In addition, the AM often takes time and considerable resources. Differences between DEM and AM results are found in predictions of the location, spacing, and dip angle of fault zones and the number of faults. Dip angles of the forward thrusts increased from the foreland to the mobile wall at the end of the experiment. Furthermore, the strain for the DEM was calculated to investigate more elaborate tectonic characteristics. These comparisons between the AM and the DEM are beneficial to the Earth Sciences community. Our results show that the DEM can successfully reproduce structures observed in the AM and also indicate the utility of the DEM in modeling largedisplacement, complex deformation of analog and geological materials. A complete structural numerical simulation laboratory with different microscopic parameters of the other materials (i.e., sand, clay, microbeads, wax, and silicone putty) can be constructed by the similar approach we presented. Our study provides a necessary comparison between commonly applied modeling approaches, e.g., discrete element simulations and scaled analog (sandbox) models, which is important for more confidently applying these methods to understand real fold and thrust belt systems.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, and further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CL and HY conceptualized the idea. CL, HY, and YZ performed the methodology. CL, CW, JZ, HY, ZW, and WW were involved in formal analysis and investigation. CL and ZW wrote the original draft. HY, ZW, DJ, SG, and RR reviewed and edited the paper. HY, CL, ZW, DJ, SG, and RR acquired the funding. (2018A-0101), and the National S&T Major Project of China (grants 2016ZX05026-002-007, 2016ZX05003-001, and 2016ZX005008-001-005). CL was supported by Program B for Outstanding Ph.D. Candidate of Nanjing University. The authors would like to thank Julia Morgan (https://earthscience.rice.edu/ directory/user/100) and Thomas Fournier for generously sharing their post-processing scripts and algorithms, which have been used to process and display the model outputs presented here. HY would also like to thank Prof. Morgan for generously sharing their discrete element code RICEBAL (v. 5.1, modified from Peter Cundall's TRUBAL v. 1.51) and Rice University for hosting his collaborative visit in 2009, providing him with the opportunity to further develop his knowledge of DEM and geomechanics principles and learn the capabilities of these methods. The authors also would like to thank Butao SHI for discussions on strain calculation and its physical significance, Chun LIU and Qian HUANG for discussions on the development of VBOX, Chuang Sun, Jian Cui, Chuan Lin, Huan LIU, and Wenqiao XU for the fruitful discussions on this manuscript, and Hui-Qun ZHOU for the maintenance of the cluster. The authors are grateful for thoughtful reviews by the editor and the anonymous reviewers, which resulted in significant improvements in the manuscript. Data used in this paper are available online at https://geovbox.com, and more examples are given in this website. Modeling results and information can be obtained by contacting CL at lichangsheng@smail.nju.edu.cn. Interested parties are encouraged to contact the authors to request copies of animated GIFs of the simulations for further study.