## 3 Model Validation  
(1) Data sources for method validation  
This research group completed an orthogonal test of oil shale-modified soil under different water content and freeze-thaw cycles. This study used the test data to verify the accuracy and error distribution of the above empirical method. For the test, a cylindrical compacted sample with a diameter of 5 cm and a height of 10 cm (dry density is 1.52 g/cm\({}^{3}\)) was used, with the same proportion and batch of the modified soil in Section ii. A half-sine pulse load was applied to the sample by the GDS dynamic triaxial tester, the load frequency was 1 Hz, and the loading duration and intermittent time of the load were 0.2 and 0.8 s, respectively. Before the formal test, the specimen was preloaded with a half-sine cyclic longitudinal load with a peak value of 66 kPa 1000 times, and a confining pressure of 30 kPa was applied. When the load was formally applied, the confining pressures of 30, 40, and 50 kPa, and the cyclic deviatoric stresses of 5, 15, and 25 kPa were used for each group of specimens as the target confining pressure and deviatoric stress amplitude of this M\({}_{\text{R}}\). The number of cyclic loads under the condition was 100. After the test, the stress and the axial recoverable deformation of the sample during the test were recorded according to the last 5 waveform curves collected under each group of stress conditions. The M\({}_{\text{R}}\) was calculated according to the following formula:  
\[M_{R}=\frac{\sigma_{0}}{\varepsilon_{0}} \tag{10}\]  
where \(M_{\text{R}}\) is the modulus of resilience, \(\sigma_{0}\) is the average axial stress amplitude (MPa) under each group of stress, and \(\varepsilon_{0}\) is the average strain amplitude (mm/mm) under each group of stress states. (2) Error distribution under different moisture contents  
The M\({}_{\text{R}}\) at different moisture contents tested in the experiment are shown in Table 3.  
In this section, the experimental data are substituted into the model (7) proposed in this study for fitting calculation, and the fitting effect and error distribution are observed. The error distribution is shown in Figure 3. \begin{table}
\begin{tabular}{c c c c c c} \hline \hline
**Confining Pressure (kPa)** & **Deviator Stress (kPa)** & **20** & **25** & **30** & **OMC (23.7)** \\ \hline  & 25 & 204.55 & 164.97 & 135.04 & 173.55 \\
50 & 15 & 228.57 & 187.77 & 165.57 & 199.57 \\  & 5 & 247.15 & 200.03 & 181.15 & 217.15 \\ \hline  & 25 & 166.77 & 120.37 & 90.71 & 136.77 \\
40 & 15 & 186.3 & 139.2 & 117.2 & 153.2 \\  & 5 & 201.78 & 157.21 & 135.78 & 174.78 \\ \hline  & 25 & 134.22 & 81.72 & 51.64 & 98.64 \\
30 & 15 & 152.22 & 104.17 & 77.17 & 120.17 \\  & 5 & 164.58 & 114.32 & 90.58 & 131.58 \\ \hline \hline \end{tabular}
\end{table} 
Table 3: Resilient modulus(kPa) of modified soil at different moisture contents. At the same time, two reference models are introduced to fit the experimental data for comparison. The first is model (1) without considering the suction contribution. The second includes the effective stress parameter proposed by [PERSON] [17] in model (1) in the same way as the model proposed in this study, which is referred to as the Khalili&K-\(\theta\) model in the following. [PERSON] [17] defines \(\chi\) as follows  
\[\chi=[u_{a}-u_{w}/(u_{a}-u_{w})_{b}]^{-0.55} \tag{11}\]  
where \((u_{a}-u_{w})_{b}\) represents the air entry value. In the fitting result, model (1) cannot fit the same set of modulus coefficients for convergence under different water contents, and the fitted coefficients for different water contents are different. The fitting statistics of model (7) and the Khalili&r-\(\theta\) model to the M\({}_{\text{R}}\) of the modified soil are compared in Table 4.  
As can be seen in Figure 3, the prediction residuals of model (7) are mainly distributed between \(-10\) and \(10\), accounting for \(85.2\%\) of the sample size, and the residual distribution is close to the normal distribution.