## 2 Nonlinear-Hysteretic Constitutive Model  
Particulate materials (i.e., rock, sandstone, sand, and sediment) exhibit strong dynamic behaviors, such as end-point memory and a closed loop (hysteresis) in the stress-strain curve, higher harmonic generation in propagating waves, a resonance frequency shift and nonlinear damping in standing waves, and slow dynamics [22; 23; 24; 25; 26; 27]. These macroscopic nonlinear responses of soils are the consequences of complex frictional mechanisms involving contact between the grains and grain rearrangement under loading-unloading cycles. The variation of the damping ratio is a physical indicator of energy dissipation as quantified by a hysteresis loop at one cycle of deformation.  
Intensive efforts have been made in the geophysics community to analytically describe these hysteretic nonlinear behaviors of granular materials. This study adopts the Davidenkov model, which has been experimentally validated, to capture the amplitude-dependent internal friction related to nonlinear hysteretic nonlinear behavior. The phenomenological equation of state was analytically derived, and its parameters were found to have physical meaning. Details can be found in the literature [28]. The equation can be rewritten in terms of the shear component:  
\[\tau=G_{\text{max}}\Big{\{}\gamma-\text{sgn}(\dot{\gamma})\frac{\alpha}{n} \Big{[}\big{(}\Delta\gamma+\text{sgn}(\dot{\gamma})\cdot\gamma\big{)}^{n}-2^{n -1}(\Delta\gamma)^{n}\Big{]}\Big{\}} \tag{1}\]  
where \(G_{\text{max}}\) is the small-strain shear modulus, \(\gamma\) is the shear strain, \(\Delta\gamma\) is the amplitude of the cyclic shear strain, \(n\) is the integer exponent related to the model order (\(n>1\)), \(\alpha\) is the nonlinear parameter, \(\dot{\gamma}(=\hat{e}\gamma/\hat{e}t)\) is the strain rate, and \(\text{sgn}(x)\) is the signum function. Numerical solutions that use a tangent formulation involve the tangent shear modulus \(G_{\text{tan}}\), which is defined as the slope of the stress-strain hysteresis loop along the loading path. By taking the derivative of Equation (1), the tangent shear stiffness is obtained:  
\[G_{\text{tan}}=\frac{d\tau}{d\gamma}=G_{\text{max}}\Big{[}1-\alpha\big{(} \Delta\gamma+\text{sgn}(\dot{\gamma})\cdot\gamma\big{)}^{n-1}\Big{]} \tag{2}\]  
It should be noted that the tangent stiffness is not associated with small strain stiffness, as the tangent stiffness \(G_{\text{tan}}\) is a mathematical concept that reveals the instantaneous state of fabric changes during a large strain test, while the small strain stiffness \(G_{\text{max}}\) is a fabric constant measurement of the stiffness [29; 30]. The backbone curve represents the trajectory of the extrema of each hysteresis curve. When the shear strain \(\gamma\) is equal to the strain amplitude \(\Delta\gamma\) in Equation (1), the piecewise function of the backbone curve is formulated for the entire strain range as follows:  
\[\tau_{\text{upper}}=G_{\text{max}}\gamma\Big{[}1-\frac{\alpha}{n}(2\gamma)^{n -1}\Big{]}\text{ for }\dot{\gamma}>0 \tag{3a}\] \[\tau_{\text{lower}}=G_{\text{max}}\gamma\Big{[}1+(-1)^{n}\frac{\alpha}{n}(2 \gamma)^{n-1}\Big{]}\text{ for }\dot{\gamma}<0 \tag{3b}\]The hysteretic nonlinear models involved with two model parameters are simulated under three different strain amplitude levels (Figure 1). The backbone curve consists of two polynomial functions, which intersect at the coordinate origin. Note that, when using the model order of \(n=2\), the first derivative of each function is equal to the small strain stiffness \(G_{\text{max}}\). In addition, higher imposed strain amplitude produces more energy dissipation and a decreasing slope (the strain-softening effect). Indeed, the hysteretic loops superimposed on the backbone curve follow the Masing rule, which has been experimentally observed in particulate materials under cyclic torsional loading [31].