## 5 Model Verification  
In order to verify the reliability of the evolution model established in this paper, the proposed model is implanted in a series of numerical studies, and the tested stress-strain curves on the GHBS samples with different hydrate saturations under various effective confining pressures in [15] are used for a comparison with the corresponding results obtained in the numerical simulations. FLAC\({}^{\text{3D}}\) software [44] developed based on the finite difference method is adopted in this numerical study. FLAC\({}^{\text{3D}}\) can be applied to solve complex problems in mechanics, especially for the mechanical responses of geomaterials. Large strain analysis can be considered in the models. Both elastic and plastic constitutive models can be implanted in the modeling, and a user-defined constitutive model is also accepted in FLAC\({}^{\text{3D}}\). Consequently, the evolution model of strength and dilation parameters can be programmed in the numerical model with FISH language. According to the experimental conditions reported in [15], numerical models of the GHBS samples are built and triaxial compression experiments are conducted under different effective confining pressures. Figure 4 presents the flowchart of the numerical modeling. A cylindrical model with 8000 meshes is generated to simulate the GHBS specimen with the size of 50 mm in diameter and 100 mm in length, as shown in Figure 5. The density values are presented in Table 7 based on [15]. The displacement boundary at the bottom of the model is fixed, and lateral stress is applied on the circumferential boundary around the model to simulate the effective confining pressure. The axial loading rate of \(-1.2\times 10^{-7}\) m/step is applied on the top boundary. \begin{table}
\begin{tabular}{c c c c c c c} \hline \hline \multirow{2}{*}{\(S_{\text{h}}\)} & \(E\) & \multirow{2}{*}{\(\
u\)} & **Initial \(c\)** & **Initial \(\varphi\)** & **Initial \(\psi\)** & **Density** \\  & **(MPa)** & \(\
u\) & **(MPa)** & **(\({}^{\circ}\))** & **(\({}^{\circ}\))** & **(\({}^{\circ}\))** & **(g/cm\({}^{3}\))** \\ \hline
0\% & 1.88 & 0.30 & 0.17 & 6.30 & 6.30 & 1660 \\
27–34\% & 2.72 & 0.35 & 0.28 & 7.65 & 7.65 & 1760 \\
41–45\% & 3.65 & 0.38 & 0.42 & 8.65 & 8.65 & 1860 \\ \hline \hline \end{tabular}
\end{table} 
Table 7: Basic parameters for the models with different hydrate saturations. Figure 4: Flowchart of numerical modeling and calculation. Figure 5: FLAC\({}^{\text{3D}}\) model of a GHBS sample. During the calculation, the Mohr-Coulomb criterion is adopted to identify whether the model is in the elastic or plastic state. The initial cohesion (\(c_{\rm i}\)) and initial internal friction angle (\(\varphi_{\rm i}\)) listed in Tables 1 and 2 are used in the criterion for the yielding condition. If the model is still in the elastic stage, the model calculation follows the elastic law with the E and \(\
u\) values presented in Table 7, as obtained from the testing results in [15]. If the model is in the yielding stage, the calculation begins to follow the plastic constitutive model with evolutionary strength and dilation parameters, as demonstrated in Sections 3 and 4, respectively. The model parameters in Tables 1-3 are used in the numerical model with the corresponding hydrate saturation. The stress-strain curves are recorded in each numerical modeling, and the curves are illustrated in Figure 6, compared with the stress-strain curves obtained in the laboratory experiments [15]. For every model with different hydrate saturations under various effective confining pressures, generally speaking, the stress-strain curves accord well with those obtained from the laboratory tests. Particularly, the strain hardening and softening behavior is simulated successfully with this model. These results exhibit the validation of the proposed model considering the evolution of strength and dilation parameters.