Introduction

The material quality indicator most commonly used for the quality of road materials is the California Bearing Ratio test (CBR). In addition, the CBR test provides a value commonly used to correlate with the resilient modulus due to its cost-effectiveness. However, this test lacks physical explanations regarding mechanical behavior because the CBR value only compares soil penetration resistance with the pressure measured in a standard material. Another issue arises from the high variability of the tests results. Then, a better understanding of the increase or decrease in the CBR value is essential, together with the variable effect of CBR value for practical design engineers or researchers in this area.

The most used relation between CBR value and the elastic modulus are linear correlations, which have great popularity, as shown by Heukelom & Klomp [1], Nielson et al. [2], and the Experimental Center for Research and Study of Building and Public Works guidelines [3]. The nonlinear correlations have gained popularity, as shown by the American Association of State Highway and Transportation Officials (AASHTO), which utilizes a nonlinear correlation ARA [4], Mallela et al. [5], Powell & Potter [6], Erlingsson [7], Leung et al. [8], Mendoza et al. [9], Gansonré et al. [10]. However, the CBR test does not directly measure its elastic modulus in coarse soils. This measure is an indirect measurement of the resistance to penetration [11]. However, it is essential to note that the literature contains a wide range of proposals with significant dispersion among them. Hence, a relationship between the CBR value and the resilient modulus is incomplete, indicating that other variables related to the behavior of coarse-grained soils may have an impact on this relationship, but they have often been overlooked.

Diverse studies show that CBR value and resilient modulus can be functions of the effects of particle shapes and sizes, drained and undrained conditions, packing structure, and compaction energy [12,13,14,15,16,17,18,19]. These effects are connected with some geotechnical parameters. For example, particle shapes and sizes can related to parameters that describe the compressibility characteristics of coarse soils [20,21,22]. The compaction energy influences shear parameters [23], hydraulic parameters (fine soils) [24], yielding stress [25]. The packing structure is related to shear behavior [26]. Nowadays, new materials are incorporated with different types of soils to improve their mechanical behavior for roads construction, such as biocatalysts, geopolymers, limes, recycled asphalt pavements, construction demolition Materials, geosynthetics, and others. Based on the research presented by Lima et al. [27], Amaludin et al. [28], Putra et al. [29], Al-Obaydi et al. [30], Li & Hao [11] has demonstrated an increment in the shear strength (i.e. undrained cohesion, cohesion, and friction angle). Although that these studies were carried out under strict laboratory conditions, there is variability between the results. One reason is the inherent variability of soils in road construction. The effects of soil variability have been addressed by diverse researchers [31,32,33,34]. However, the variability of geotechnical properties in road construction and its importance have only been examined for the asphalt mixtures, deflections in pavement structures, spatial variation of laboratory tests, and quality of materials, e.g. [35,36,37,38,39,40,41]. Then, the variability effect needs to be understood in complete form.

To address the challenges mentioned above, this paper addresses CBR behavior in coarse-grained soils. It uses Finite Element Method (FEM) simulations in which geotechnical parameters vary randomly. These parameters may be related to the previously mentioned effects. Two constitutive models are used. One is an elastoplastic model with an elastic part and a Drucker-Prager rupture criterion with a cap. The other has an elastic part and a Mohr–Coulomb rupture criterion. The capped model takes into account the compressibility characteristics of the soil. Monte Carlo analyses were performed on the geotechnical parameters using both models for the CBR test. Finally, a sensitivity analysis of each parameter and state variable was conducted to assess their influence on predicting the CBR value.

Based on the analyzes conducted, it was shown that several parameters are crucial for understanding the CBR value, particularly when taking into account the variability effect through the penetration-stress curve. The increase or decrease in the CBR value is a function of the elasticity modulus, yielding stress and friction angle. Besides, the simulations expand the knowledge of the shearing mechanisms, generated stresses, displacement fields, and load sharing when the CBR test is carried out. From these results, a physical explanation of test results was obtained. FEM simulations showed stress paths (or zones) in elastic, compression, and shear behavior, explaning the importance of elastic modulus, yielding stress, and friction angle in the CBR results. In addition, with analyzes conducted, a new equation was proposed and compared with practical equations proposed by international standards, recommendations, and other sources to determine the probability of underestimated values CBR according to the correlations used. Overall, the analysis provides valuable insights into the variables (or parameters) that result in either high or low CBR values for coarse-grained soils and which practical equations can provide more accurate results.

Literature review on factors and modeling that influence CBR value for coarse-grained soils

The California Bearing Ratio (CBR) test is a straightforward procedure that involves penetrating a soil sample using a standard 2-inch (5.08 cm) punch at a controlled speed of 1.25 mm/min. The test measures the soil's response to penetration, and a stress-penetration curve is obtained up to a depth of 12.5 mm. The stresses recorded at 2.5 and 5 mm penetrations are normalized to the standard stress values obtained on the reference material. The standard stress used for normalization is 6900 kPa at 2.5 mm penetration and 10,300 kPa at 5.0 mm penetration. The reference material historically used in the CBR test was limestone, as it demonstrated favorable results during the test's development in the 1920s in California. However, many other materials worldwide now exhibit better characteristics in terms of strength and performance. Given its simplicity, the CBR test has become the most commonly used method globally for deciding the thicknesses of pavement structure layers and the quality of bases, sub-bases, and subgrades. Finally, the primary objective of the test was to provide an efficient means of measuring soil strength under the wheels of heavy vehicles, which is why it has gained popularity in the field.

In geotechnics of road, the CBR test is typically approached as an elastic problem, considering the elastic modulus. Magnan and Ndiaye [18] assumed a conical stress distribution below the punch. They further considered a homogeneous stress distribution on a cylinder beneath the cone. As a result, the total elastic displacement of the punch is obtained by summing the strains in the conical and cylindrical regions. Based on this approach, the following Young's modulus was derived [42]:

$$E=\frac{{p}_{m}d}{{h}_{e}D}\left(H+\frac{d\left(L-H\right)}{D}\right)$$
(1)

where E is young’s modulus, he is the elastic displacement of the punch, pm is the mean stress under the punch, d is the diameter of the punch, D is the diameter of the CBR mold, L is the height of the mold of the specimen tested, and H is the height of the cone. Magnan and Ndiaye [18] assumed that a 45◦ angle is close to the friction angle of the material, and used it as the opening angle of the cone. When the stress and indentation corresponding to CBR = 100% (pm = 6.9 MPa and Δhe = 0.254 cm) are used together with geometry of the mold given by: d = 4.94 cm, D = 15.24 cm, L = 12.7 cm, and H = 5.15 cm, Eq. 1 results in E = 96.3 MPa. So, the relationship between the CBR and Young’s modulus is linear. The analysis based on contact theory is only valid for the elastic case. On the other hand, Nielson et al. [2] proposed a correlation for CBR based on classical soil mechanics principles, as follows:

$$E(kPa)=\frac{0.75\pi a{\left(1-\nu \right)}^{2}}{\left(1-2\nu \right)}*689.5CBR$$
(2)

Assuming \(\nu\) (Poisson´s ratio) = 0.25 and a (indenter radius) = 2.48 cm (or 0.975 inches), the equation changes to E(kPa) = 1800CBR. However, based on experimental evidence, Nielson et al. [2] modified the equation’s constant from 1800 to 2150 [43]. On the other hand, some empirical and semi-empirical formulas have also been proposed to estimate the CBR value or assess the behavior of road materials. Heukelom and Klomp [1] proposed using E(kPa) = 10340CBR; the guidelines of the French Experimental Center for Research of Building and Public Works [3] proposed E(kPa) = 5000CBR, and the American Association of State Highway and Transportation Officials (AASHTO) proposed a nonlinear equation: E(MPa) = 17.6CBR0.64 [4, 5]. Mendoza and Caicedo [19] proposed a nonlinear correlation between CBR and E with the inclusion of pre-consolidation stress (pP) and compression index (\(\lambda\)), as shown below:

$$CBR[\% ] = - 1.3 \times 10^{{ - 9}} \lambda ^{{0.236}} E^{2} + \left( {0.002\;ln\;p_{p} - 0.00086} \right)E$$
(3)

The purpose of the present study is to analyze the process of indenting soil with a cylindrical punch, as in the CBR test by FEM models. These FEM models were developed using random parameters associated with coarse-grained soil.

On the other hand, the effects of particle shapes and sizes in the CBR test were investigated by Zhang et al. [26], who made some micromechanical simulations to change the particle shapes and sizes in the biaxial and CBR tests. Similar experimental results were obtained by Tutumluer et al. [14], Mishra et al. [15], and Kwon et al. [16], where different types of material were tested with different types of particle shapes, sizes, and breaking. The gradation change in CBR value with different compaction energies was studied by Magnan and Ndiaye [18]. Other research on energy compaction was conducted by Benson et al. [24] and Sreelekshmypillai and Vinod [17] in fine soils. Araya et al. [23] showed the influence of compaction energy in shear stresses in CBR and triaxial tests. Then, all these effects are related to geotechnical parameters. Mesri and Vardhanabhuti [44] show in one-dimensional compression that the compression index and yielding stress (the point of maximum curvature) are influenced by the shape and composition of the material grains. Clements [45] and Oldecop and Alonso [46] demonstrated that angular particles are more prone to breakage than rounded particles when subjected to equal loads. Other research report abrupt slope changes during one-dimensional compression due to particle fracture, splitting, and rearrangement [20,21,22]. Note that the yielding stress and compression index are parameters in one of the models studied later. Araya et al. [23] conducted triaxial tests with compaction energy increments, showing increased cohesion and friction angle. Additionally, some unconventional materials are mixed with soil to improve the CBR value and geotechnical parameters. Tamassoki et al. [47] mixed residual soil, coir fiber, and lime to improve the CBR value, cohesion, elastic modulus, and friction angle with positive results. Similar results were found by Al-Obaydi et al. [30] with mixed of the soil and construction-demolition material. Lima et al. [27] address the mix of recycled asphalt pavement, sedimentary soil from Brazil, and Portland cement. These mix tests show improvement in the CBR value and unconfined compressive strength (correlated with undrained cohesion).

The incorporation of parameter or state variable variability has been a focal point in numerous investigations conducted [48,49,50,51,52]. Some recent studies worth mentioning in the present context of transportation geotechnics are the following: Caro et al. [35], and Castillo et al. [36] have studied the effect of variability in the deterioration of asphalt mixtures. Thyagarajan et al. [37] studied the vertical surface deflections in pavement structures. They aimed to establish a relationship between deflection-basin-related indices and load-induced responses in pavement structures while considering the variability effect of surface deflections. The researchers employed various techniques such as normal distributions, lognormal distributions, Monte Carlo analysis, random fields, among other methods. In relation to the CBR value, a study by Look [38] takes into account the spatial variation effect in order to determine the suitable design CBR value. The results show that the lognormal distribution is an appropriate statistical model to characterize the design CBR. Erzin and Turkoz [39] used neural networks to predict the CBR value for granular soils based on the input of specific gravity, coefficient of uniformity, coefficient of curvature, dry density, water content, and mineralogic contents. Based on these analyses, a predictive model for the CBR value was developed. Trong et al. [40] conducted a sensitivity analysis using three soft computing models to predict the CBR value in soils containing fines. The sensitivity analysis used 10 input parameters: gravel percentage, coarse sand percentage, fine sand percentage, fine soil material, organic matter content, liquid limit, plastic limit, plasticity index, optimum moisture content, and maximum dry density. This analysis was based on 214 soil samples. Also, Kurnaz and Kaya [41] conducted analyses using neural networks and multiple linear regression models for soils in Turkey. They developed several prediction models to determine the CBR value of compacted soils based on their analyses.

Numerical simulation of a CBR test

Constitutive models used in simulations of the CBR test and materials

Throughout the XVII century, several yield criteria were developed within a mathematical framework for different materials. These criteria have demonstrated excellent adaptability to different materials. However, not all materials or conditions can be simulated using the same criteria. For example, the Tresca and von Mises criteria do not account for the effects of isotropic stress (i.e., the mean stress), while the Mohr–Coulomb and Ducker Prager criteria consider the increase in strength produced by the mean stress [53,54,55]. In this study, finite element simulations were conducted using an elastic model that incorporates multiple yield criteria. The elastic part of the model depicts the material's behavior up to the point where the yield criterion is reached, and elastoplastic strains begin to occur. The relationship between stress and strain in the elastic region is described by an elastic tensor that depends on the elastic modulus E and Poisson's ratio ν. Upon reaching the yield criterion, plastic strains are introduced as the stress paths intersect the criterion's limit. In this work, two yield criteria were employed. The first criterion was the Mohr–Coulomb, whereby the behavior of stresses paths is determined based on the following equation:

$$\tau =c+{\sigma }_{n}\text{tan}\varphi$$
(4)

where \(\tau\) is the shear stress of the soil, c is the cohesion of the soil, \({\sigma }_{n}\) is the normal stress, and \(\varphi\) is the friction angle of the soil. Equation 4 shows the shear strength increase as normal stress increases due to cohesion (or the influence of isotropic stress), as shown in Fig. 1a. One advantage of this criterion is that parameters are common in transportation geotechnics. Therefore, the elastoplastic model only needs four parameters to characterize soil behavior accurately. Additionally, the dilatancy effect, which is significant in granular materials, can be described using a non-associated flow rule. This importance of this dilatancy effect was demonstrated by Strahler et al. [56] and Bolton [57]. To describe the dilatancy effect, the dilation angle was simplified using an empirical relationship proposed by Bolton [20] as ψ = φ—33°. For φ > 33°, the dilation angle is equal to zero. The dilation controls the plastic volumetric strain developed during plastic shearing. This dilatation angle should be modeled as a function of the soil strain. However, for the sake of simplicity, a constant value was employed in this study.

Fig. 1
figure 1

Scheme of the shear criteria: a Mohr–Coulomb criterion, b Drucker-Prager criterion with cap (adapted from [9, 55])

Full size image

The second criterion employed in this study is the Drucker-Prager with an elliptical cap (Fig. 1b). The equation of this criterion is similar to the previous criterion:

$$t = d + p\;{\text{tan}}\beta$$
(5)

However, the equation factor changes to other similar factors, as demonstrated below. The cohesion factor is denoted as d, the mean effective stress is represented by p, and the friction angle in a p–t plane (as shown in Fig. 1b) is symbolized as\(\beta\). Furthermore, the shear stress, measured by t is defined in the following equation:

$$t=\frac{q}{2}\left[1+\frac{1}{K}-\left(1-\frac{1}{K}\right){\left(\frac{r}{q}\right)}^{3}\right]$$
(6)

where q is the deviator stress, r is the third stress invariant, and K is a constant that governs the shape of the yield surface in the plane \(\pi\). Equation 6 can be written as follows:

$$t=\frac{q}{g}$$
(7)

with,

$$g=\frac{2K}{1+K+\left(1-K\right){\left(\frac{r}{q}\right)}^{3}}$$
(8)

Parameters d and β in the Drucker-Prager criterion can be related to the friction angle \(\phi\) of the soil, and the cohesion c for the Mohr–Coulomb criterion using the following equations:

$${\text{tan}}\beta = \frac{{6\;{\text{sin}}\varphi }}{{3 - {\text{sin}}\varphi }}$$
(9)
$$d = \frac{{18c\;{\text{cos}}\varphi }}{{3 - {\text{sin}}\varphi }}$$
(10)

However, if K is less than one, the yield stress in tension and compression in triaxial conditions are not the same. The range to k is 0.778 to 1 to ensure that the yield surface remains convex.

The elastoplastic model based on the Drucker-Prager criterion and a cap ha incorporates s a cap-yielding surface characterized by three associated parameters. These parameters include the eccentricity of the cap ellipse (Fig. 2), and a constant that governs the shape of the yield surface in the plane \(\pi\) (the plane orthogonal to the hydrostatic pressure axis in the principal stress space). Another important variable is pp which serves as a state variable that describes the evolution of the material's isotropic compression law (which controls hardening softening). This stress corresponds to the yield point on the isotropic compressibility line. The isotropic compression law is influenced by volumetric plastic strains \(\varepsilon_{v}^{p}\) where λ is the slope of the virgin isotropic compressibility line in the plane constituted by the void ratio e and the natural logarithm of p; \(\kappa\) is the slope of the unloading–reloading line of the material, e0 is the initial void ratio, as shown below.

Fig. 2
figure 2

Scheme of the yield surface of the Drucker and Prager criterion: a on the p-q plane, b on the π-plane (adapted from [57, 58])

Full size image
$${\varepsilon }_{v}^{p}=\frac{\lambda -\kappa }{1+{e}_{0}}\text{ln}\frac{p}{{p}_{p}}$$
(11)

The previous model, which incorporates shear strength and compressibility, can be described using nine parameters. These parameters can be determined using classical geotechnical laboratory tests by means of triaxial tests and compression tests.

The simulations used coarse-grained soil that showed drained behavior due to its high hydraulic conductivity. The granular material selected for the simulations was crushed hornfels rock as reported by Araya et al. [23] sourced from a quarry in South Africa. This material was chosen due to its comprehensive characterization the availability of all the necessary test data from the same study. Other materials reported in CBR tests include uncrushed gravel, limestone, and river sand described by [14, 15, 59]. However, it is worth noting that the characterization of these materials was carried out with specific purposes in mind, which may differ from the requirements of the current study or simulation. Table 1 presents the initial parameters and variables used for crushed hornfels rock. The typical coefficients of variations of geotechnical parameters shown in Table 1 were found in Kirby [60], Phoon and Kulhawy [31], Griffiths et al. [32], Llano-Serna et al. [33] and Zevgolis et al. [34]. Note that in this study, it is assumed that the parameters are not correlated with other parameters or relations. One of the objectives of this study was to observe the individual influence of each parameter on the CBR value. The parameters are defined by relationships (i.e., models) that describe specific physical situations using constants, following the approach of Cividini et al. [61] and Lei et al. [62]. These constants (or parameters) often represent the inherent properties of materials. While there may be statistical relationships between parameters, this paper did not consider them. Some recently published papers that address the dependent parameters include works by Brinkgreve [63] and Bolaños and Hurtado [64].

Table 1 Initial parameters and variables used
Full size table

Finite element analyses

Finite element (FEM) analyses of the CBR test were conducted to assess the significance of variability in geotechnical parameters during the test and to evaluate the suitability of the rheological models in accurately simulating the physical phenomena. The initial steps of the simulation and the final calibration provided valuable insights into the underlying shearing mechanisms and the generation of stresses.

The Abaqus environment was used to enable the 3D analyses of the CBR test. The standard CBR test mold is 6 inches (15.24 cm) in diameter and 7 inches (17.78 cm) high with a 2-inch (5.08 cm) ring. The same geometry was used in the FEM models (cylindrical sample with 15.2 cm diameter and 11.65 cm height). Lateral confinement of the soil was achieved by using a CBR mold wall, and the interaction between the soil the wall was modeled using friction interaction. Sensitization analyses were conducted by varying the friction factors, as shown in Fig. 3a. The friction coefficient was calculated as a percentage of the friction angle (\(\mu\) = ĉtan\(\phi\)) of the soil. The simulations were conducted using the following values for the coefficient of friction ĉ = 1, 0.875, 0.75, 0.625, 0.5, and 0.325. The initial friction angle was set to 23°. In the model base, a fixed boundary condition was applied. Figure 3b provides a schematic representation of the CBR geometry. Note that the initial simulations assumed a low value of soil permeability. The following analyses (Sect. "Simulations result with random parameters in FEM models") used k = 0.1 m/s to simulate a drained condition.

Fig. 3
figure 3

Finite Element model (FEM) for the CBR test (from Abaqus 6.21)

Full size image

To observe response stability and simulation time, the number of elements in the model was adjusted based on the described geometry, as illustrated in Fig. 4a. The vertical stress and pore pressure were extracted at the center point of the model on the surface, as shown in Fig. 4b. Note that this point remains consistent across all models. Stresses and pore pressure were obtained at a strain of 12.5 mm. Additionally, different types of elements were tested in the FEM model to investigate the resulting stress and strain fields, as well as the simulation time. In this study, the soil was modeled using C3D8P (continuous, eight nodes with pore pressure measurement), while the mold wall was represented by C3D8R elements (continuous, eight nodes with reduced integration).

Fig. 4
figure 4

a Simulation time as a function of the number of elements, b Stress and pore pressure response at a specific point for different numbers of elements

Full size image

In the model, 12,096 elements were utilized, as depicted in Fig. 4b. Subsequently, a modification was made to the bias ratio. The bias ratio refers to the difference in length between nodes at opposite ends of an edge. This adjustment was implemented to position smaller elements near the test piston and larger elements near the mold walls. Figure 5 demonstrates that a higher bias ratio yields a response similar to that achieved by increasing the number of elements.

Fig. 5
figure 5

Change of response in the center point with the change of bias ratio

Full size image

The CBR test consists of compacting a material to a certain energy and water content (obtained from the Modified Proctor Test) and then immersing the soil-filled mold for four days with a 4.5 kg overweight on top to simulate pavement layers. The soil sample in the mold is then penetrated by a standard 2-inch (5.08 cm) punch at a speed of 1.25 mm/min. Simulations were performed in two sequential stages within the models. First, the pavement construction stage involved placing a 4.5 kg overweight on top of the FEM model. Next, a vertical boundary condition was imposed to induce strain in the plunger area of the model (simulation of the cylindrical punch used to CBR test). The boundary condition was set to maintain a constant strain rate of 1.25 mm/min until a total displacement of 12.5 mm was achieved. This boundary condition assumes that the cylindrical punch is much more rigid than the soil.

Monte Carlo analysis for CBR value

With the constitutive models and FEM models assumed to be adequate, the CBR value can be prognosticated using Monte Carlo analysis. The analysis proceeds in the following steps: the distribution functions of the input parameters are studied, and means and standard deviations are estimated. In general, the input parameters are assumed to be independent, and samples are taken from the marginal distribution of each factor [65]. However, parameters may exhibit statistical relationships, which was not considered in this paper. Some publications discuss the dependent parameters [63, 64]. For each of the input parameters (see Table 1), a sample is extracted, producing a set of row vectors. These samples are then combined to form a matrix, representing the total number of simulations

$${\varvec{M}}=\left[\begin{array}{ccccccc}{E}^{\left(1\right)}& {c}^{\left(1\right)}& {\phi }^{\left(1\right)}& {\nu }^{\left(1\right)}& {\lambda }^{\left(1\right)}& {\kappa }^{\left(1\right)}& {{p}_{p}}^{\left(1\right)}\\ {E}^{\left(2\right)}& {c}^{\left(2\right)}& {\varphi }^{\left(2\right)}& {\nu }^{\left(2\right)}& {\lambda }^{\left(2\right)}& {\kappa }^{\left(2\right)}& {{p}_{p}}^{\left(2\right)}\\ {E}^{\left(3\right)}& {c}^{\left(3\right)}& {\varphi }^{\left(3\right)}& {\nu }^{\left(3\right)}& {\lambda }^{\left(3\right)}& {\kappa }^{\left(3\right)}& {{p}_{p}}^{\left(3\right)}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ {E}^{\left(N-1\right)}& {c}^{\left(N-1\right)}& {\varphi }^{\left(N-1\right)}& {\nu }^{\left(N-1\right)}& {\lambda }^{\left(N-1\right)}& {\kappa }^{\left(N-1\right)}& {{p}_{p}}^{\left(N-1\right)}\\ {E}^{\left(N\right)}& {c}^{\left(N\right)}& {\varphi }^{\left(N\right)}& {\nu }^{\left(N\right)}& {\lambda }^{\left(N\right)}& {\kappa }^{\left(N\right)}& {{p}_{p}}^{\left(N\right)}\end{array}\right]$$
(12)

N is the simulation number from the Monte Carlo experiment. This number can be obtained when the objective value's mean and standard deviation are stabilized (CBR value), as shown by Bolaños and Hurtado [64], Haldar and Babu [66], Al-Bittar et al. [67], among others. The vector C can be obtained by computing from the FEM model with the matrix M. This vector contains the desired CBR values for N simulations, as illustrated below.

$${\varvec{C}}=\left[\begin{array}{c}{CBR}^{\left(1\right)}\\ {CBR}^{\left(2\right)}\\ {CBR}^{\left(3\right)}\\ \cdots \\ {CBR}^{\left(N-1\right)}\\ {CBR}^{\left(N\right)}\end{array}\right]$$
(13)

The vector C can be utilized to conduct a sensitivity analysis, which provides valuable insights into the uncertainty distribution and its impact on the CBR value. Various tools can be employed in this analysis, including scatterplots and slicing of the scatterplots. These tools assist in obtaining sensitivity measures and identifying the main parameters influencing the CBR value, as shown below.

Generation of random numbers and integration into the CBR model

The purpose of this paper is to shed light on the influence of random geotechnical parameters on the CBR value. To gain a comprehensive understanding, random numbers were generated for the geotechnical parameters listed in Table 1. The random numbers were obtained by following a log-normal distribution. Monte Carlo simulations were then conducted using the FEM model (see Sect. "Monte Carlo analysis for CBR value"), incorporating these geotechnical parameters. The results highlight the importance of geotechnical parameters in influencing the stresses-penetration curves, thereby contributing to the physical simulation of the CBR test.

Random finite element analyses were performed as follows:

  1. a.

    1050 random numbers were generated for each geotechnical parameter of the elastoplastic model (with Drucker-Prager failure criteria), and a simulation was performed for each number (1050 simulations in total). This number was used because it stabilized the mean and the standard deviation of the objective value (CBR value) [64, 65, 67]. The random numbers were generated using statistical variables as shown in Table 1 for a lognormal distribution [31, 68].

  2. b.

    The parameters of the finite element model were set by using a subroutine in python language and processing it automatically. This subroutine is performed with a nested loop for shipping multiple processes simultaneously in order to obtain better performance in computational terms;

  3. c.

    The important output variables of the problem (strains and stresses) were saved in an external file.

  4. d.

    The results were analyzed until a penetration of 12.5 mm was reached, and all parameters obtained from the simulations were saved.

  5. e.

    The CBR value at a penetration of 0.0025 m or 0.005 m was obtained from the stress-penetration curves..

Simulations result with random parameters in FEM models

The results of the CBR simulations were plotted using random parameters, as depicted in Fig. 6, for different criteria. Figure 6a shows the stress-penetration curves for the Drucker-Prager criterion. The results are presented for stress-penetration curves at a penetration of 0.0025 m. Figure 6b displays the histograms of the CBR values, along with the calculated lognormal density curve. The cumulative probabilities indicate that the CBR simulations follow a lognormal distribution, as validated by a Kolmogorov–Smirnov test with a significance level of 5%. The standard deviations and mean values for different penetration depths (0.0025 and 0.005 m) and criteria are shown in Table 3.

Fig. 6
figure 6

a Comparisons between simulations and penetration to 0.0025 m., b Density curve and histogram of CBR for simulations for Mohr–Coulomb and Drucker-Prager criterion

Full size image

The results in Table 2 show an increase in CBR value with increasing limit penetration for both models, with respect to a mean value and a standard deviation. However, in all four cases, the coefficient of variation remains almost constant. The increase in the mean value with the Mohr–Coulomb criterion was 12.75% compared to the model with the Drucker-Prager model at 0.0025 and 31.75% at 0.005 mm. The comparison indicates that the increase in the CBR mean value with the penetration increase was 9.40% for the Drucker-Prager model and 27.84% for the Mohr–Coulomb model. The standard deviation variation was 6.40% for the Drucker-Prager model and 32.04% for the Mohr–Coulomb model.

Table 2 Statistical variables from simulations for different penetrations for CBR values
Full size table

Capacity of elastoplastic models for representation CBR for different materials

FEM simulations were utilized to depict the behavior of coarse-grained soil during the CBR test. These simulations incorporated variability effects in the geotechnical parameters to accurately represent crushed hornfels rock, as reported by Araya et al. [23]. The initial step involved an extensive literature review focused on granular materials with a complete characterization from a mechanical behavior standpoint. It became evident that while the literature contains many reported materials, they often lack complete characterizations due to the specific intended use of the tests conducted. The selected materials were crushed hornfels rock (CHR) [23]; uncrushed gravel (UCG) [14, 15]; limestone (LM) [14, 15]; and river sand (RS) [59]. These materials are used to construct roads in South Africa, the United States, and Korea. Figure 6 shows the capacity of the two elastoplastic models selected to represent the CBR value for different materials. However, the random parameters were generated based on the mean parameters of crushed hornfels rock, as reported by Araya et al. [23]. Furthermore, Fig. 6 was generated based on the elastic modulus, as it is the parameter most commonly used to establish correlations with the CBR value.

Drained and undrained response

Figure 7 shows the pore pressure generation in the models with different permeabilities. Figure 7a shows a low pore pressure in the FEM model. This model was made with a k = 0.1 m/s. Figure 7b shows a high pore pressure with k = 1.80E-09 m/s. This is important because the CBR test does not control the drain's initial condition (drained and undrained condition). Thus, coarse-grained soil has an initial drained condition, and finer soil may present an initial undrained condition. This condition is a function of the permeability value of the material and the compaction energy. In addition, the FEM models can represent drained and undrained conditions generated in the CBR test. The simulations presented in Figs. 3, 4, and 5 were conducted using a low permeability value. However, for the remaining figures, a high permeability value (k = 0.1 m/s) was employed to better represent real conditions typically encountered in coarse materials.

Fig. 7
figure 7

Pore pressure generated in the Finite Element model for: a drained condition; b undrained condition

Full size image

Simulation comparisons of elastoplastic criteria

Figure 8 illustrates the distribution of stresses in the coarse-grained soil for two criteria considering random parameters. The results indicate the presence of three distinct stress zones. The first zone corresponds to high-shear stress, where the shear criteria come into play. In this zone, the rupture criteria operate based on the shear parameters considered in the two criteria (see Sect. "Constitutive models used in simulations of the CBR test and materials"). The second zone had high compressive stresses. The Drucker-Prager criterion with a cap can mobilize the surface of the elastic domain. This mobilization is a function of yielding stress in the compression paths. The effect of the isotropic stress cannot be represented with an elastoplastic model with the Mohr–Coulomb criterion without a cap. A cap is a real effect observed in geomaterial, as is discussed in the next section. The third zone corresponds to an elastic zone, which is defined by Young’s modulus and Poisson’s ratio. Mendoza and Caicedo [19] compared the stress paths using an elastic model and an elastoplastic model in different zones of a FEM model for CBR. The stress paths showed that many paths stay in elastic states; some reach the cap, and only a few reaches the shear criterion. Figure 8b shows high-concentration stresses reached with the Mohr–Coulomb criterion compared to Drucker–Prager. This is due to the fact that the simulations performed with the Mohr–Coulomb do not consider a cap. So, the stress paths below the shear criterion are not limited by yielding stress. This behavior does not change with the influence of variability in geotechnical parameters.

Fig. 8
figure 8

FEM analysis contour and surface at stresses in the penetration of 2.5 mm

Full size image

The CBR test shares similarities with indentation tests employed for characterizing metals. In the case of a cylindrical punch penetrating an elastic half-space, it represents a non-Hertzian contact problem. This scenario results in a significant stress concentration at the edges of the cylindrical punch. Sneddon [69] derived the solution for the pressure distribution underneath the punch during penetration. This solution is only valid for an elastic case. Numerous attempts have been made by researchers to obtain analytical elastoplastic solutions for indentation problems, yet none of these proposals have successfully captured the stress field generated by a cylindrical punch, as demonstrated in the works of Johnson [70] and Riccardi & Montanari [71]. Despite these limitations, the elastoplastic behavior offers notable advantages as it allows for the consideration of yield stress, hardening mechanisms, and other soil characteristics. Consequently, several studies have employed Finite Element Method (FEM) simulations to represent indentation with elastoplastic behavior [19, 72,73,74]. However, it is important to mention that none of these studies explicitly address the effect of variability in their analyses.

Sensitivity analysis of the input parameters predicting the CBR value

A sensitivity analysis was conducted using FEM simulations with elastoplastic criteria to examine the relationship between geotechnical parameters (inputs) and CBR values (outputs). The aim was to identify the most important geotechnical parameters in the CBR value using various methods. These analyses were carried out according to Mohr–Coulomb and Drucker-Prager criteria. In addition, the CBR value was obtained for penetrations of 2.5 and 5 mm.

The first method employed was Scatterplots versus Derivatives. Figure 9 illustrates the scatterplots of geotechnical parameters (GPi) against CBR values obtained through Monte Carlo simulations (refer to Sect. "Monte Carlo analysis for CBR value"). Scatterplots are a type of data visualization that depicts the relationship between two numerical variables. Through these scatterplots, distinct patterns were observed, particularly in relation to Young's modulus and yielding stress for the Drucker-Prager criterion, as well as Young's modulus and friction angles for the Mohr–Coulomb criterion. However, the straightforward derivative (SS = ∂CBR/∂GPi) cannot be used in the present analyses because the geotechnical parameters have different orders of magnitude. Thus, the parameters were normalized by dividing each parameter by its respective maximum value. By doing so, the normalized derivatives (SSN = ∂CBR/∂(GP/GPmax)) were obtained as shown in Table 3, where a higher SSN value indicates a stronger pattern. Utilizing Scatterplots versus Derivatives proved to be a straightforward and informative approach for conducting the sensitivity analysis [65].

Fig. 9
figure 9

Capacity of elastoplastic models to represent CBR test

Full size image
Table 3 Constant of the fit linear equation and sensibilization parameter for all the cases studied in CBR value
Full size table

The second method employed the use of the response vector C (output vector) and matrix M containing geotechnical parameters. It is important to highlight that the C output vector's size corresponds to the Monte Carlo experiment (N = the number of times that FEM simulations run). An example of this is shown in Fig. 9. Subsequently, a straightforward linear regression analysis is applied to the collected data, as shown in Fig. 9 and Table 3. This method has gained significant popularity [65]. The equation for linear regression is shown below:

$$CBR\left(i\right)={b}_{0}+\sum_{j=1}^{r}{b}_{zj}{z}_{j}^{i}$$
(14)

The coefficients \({b}_{0}\) and \({b}_{zj}\) are determined through a least-square computation, as explained in the paper. It is crucial to note that the sample size is sufficiently large to ensure the representativeness of these values (see Sect. "Generation of random numbers and integration into the CBR model"). In addition, the coefficient \({b}_{0}\) is represents the intercept in the linear equation, \({b}_{zj}\) is the slope, and \({z}_{j}^{i}\) denotes a specific geotechnical parameter for a given simulation. The results for these coefficients are shown in Table 3. Table 3 displays the results for these coefficients, including the value of R2. In this context R2 represents the square of the sample correlation coefficient (r) between the observed outcomes and the observed predictor values for simple linear regression [75]. The computed R2 value is used to assess the prediction of future results or the testing of hypotheses based on other related information. It represents the proportion of the variation in the dependent variable that can be predicted from the independent variable(s). R2 ranges between 0 and 1 where a value close to 1 indicates a strong relationship between the variables and a high predictability of the dependent variable. However, this assumption is only valid for linear regression regression models, as they only consider the linearity of the observed relationship. Therefore, the degree of linearity can decompose the variance of CBR, even in cases where the R2 value is low, allowing for sensitive analyses.

The third method is sigma-normalized Derivatives, where the derivative is normalized by the input–output standard deviations, as shown below:

$${S}_{{GP}_{i}}^{\sigma }=\frac{{\sigma }_{{GP}_{i}}}{{\sigma }_{CBR}}\frac{\partial CBR}{\partial {GP}_{i}}$$
(15)

where \({\sigma }_{{GP}_{i}}\) and \({\sigma }_{CBR}\) are standard deviations of geotechnical parameter (input) analyses and CBR (output). Furthermore, it is important to note that the sum of squared terms in the analysis is equal to one \(\sum {\left({S}_{{GP}_{i}}^{\sigma }\right)}^{2}=1\). This sensitivity measure \({\left({S}_{{GP}_{i}}^{\sigma }\right)}^{2}\) is consistent with Fig. 9, shown in Table 3. This measure is more convincing than the other two methods because the relative ordering of the geotechnical parameters now depends on both slope (first method) and the relationship between the standard deviations of input–output. Also, the sensitivity measures are neatly normalized to one. However, it is important to note that in order to obtain a value equal to one, the number of simulations must be sufficiently large. This recommendation for sensitivity analysis comes from Intergovernmental agencies [65].

Sensitivity analysis serves the purpose of preventing errors that may arise from mistakenly designating a non-influential factor as important or incorrectly categorizing an essential factor as non-influential. Saltelli et al. [65] conducted a comprehensive discussion addressing these errors, along with other potential types of errors.

According to Table 3, the parameters that have the most influence on the CBR (California Bearing Ratio) value in the elastoplastic model with the Mohr–Coulomb criterion are Young's modulus and friction angle. On the other hand, in the elastoplastic model with the Drucker-Prager criterion with a cap, the most important parameters affecting the CBR value are Young's modulus, isotropic yield stress, and friction angle, as indicated in Fig. 9. These results were based on three sensitivity parameters (SSN, R2 and \({\left({S}_{{GP}_{i}}^{\sigma }\right)}^{2}\)), where a high value of the parameter corresponds to a greater influence on the dependent variable. The CBR value depends not only on Young's modulus, as typically used in correlations, but also on the soil's isotropic yielding stress and friction angle, making it important to consider these factors for accurate assessments.

The findings presented in Table 3 are consistent with the numerical results reported by Caicedo and Mendoza [76], Shaban and Cosentino [77], and Narzary and Ahamad [74], which also employed the Mohr–Coulomb criterion without a cap. However, the results from Mendoza and Caicedo [19] demonstrate that the presence of a cap has a significant impact on the behavior of the CBR value. This can be explained by the fact that stress paths predominantly remain within the elastic domain, with some stress paths mobilizing to the cap yield surface when using the Drucker-Prager criterion. As a result, the influence of parameters in the model shifts, with the CBR value calculated with the Drucker-Prager (with a cap) showing an increased effect of yield stress and friction angle, while the influence of Young's modulus decreases. Consequently, the results shown in Table 3 and Fig. 9 exhibit similar trends, although the influence of parameters differs when calculating the CBR value with a 5 mm penetration using the Drucker-Prager criterion.

Discussion of results and comparison between elastoplastic simulations and correlations for practical design

The most important parameters in the simulations were Young’s modulus, yielding stress, and friction angle. As expected, Young’s modulus was the most important parameter in the behavior of CBR value. However, higher-yielding stress causes an increase in the CBR value. This increment is not linear as yielding stress increases. The yielding parameter in granular materials is strongly related to crushability which, in turn, depends on the shapes of particles, their strength, and the grain size distribution [13,14,15,16]. In the Drucker–Prager criterion with a cap, the friction angle exhibits a relatively small effect compared to the Mohr–Coulomb criterion. However, the two models studied show an effect of the friction angle. Some experimental and numerical research showed that the friction angle has a minor impact on compacted materials with high friction angles. For materials with low friction angles and high moduli, there are significant reductions in the CBR value [19, 23]. Based on the numerical simulations of this study, the stiffer materials reach the failure criterion at an earlier stage, specifically at a penetration of 2.5 mm, which is the point at which the CBR value is calculated. Nonetheless, it is worth noting that an extremely high friction angle (close to 80°) is unrealistic for most materials.

As was shown in Sect. "Literature review on factors and modeling that influence CBR value for coarse-grained soils", yielding stress exhibits various effects, such as changes in slope during one-dimensional compression, particle fracturing, splitting, and particle rearrangement, and the shape and composition of the material grains [20,21,22, 44,45,46]. Note that the yielding stress observed in one-dimensional compression is analogous to the yielding stress (pP) in isotropic compression by the earth pressure coefficient at rest, representing the size of the boundary surface in the Drucker-Prager model with a cap (see Fig. 2a). Experimental evidence supports the influence of yielding stress on the yield surface, with the size of the surface directly correlating to the magnitude of the yielding stress. The influence of yielding stress on the yield surface has been observed in various geomaterials, as demonstrated by studies such as Roscoe et al. [78], Tavenas et al. [79], Graham et al. [80], Rutter and Glover [81], Alonso et al. [82], among others. The Drucker–Prager model with a cap incorporates the influence of yielding stress within the cap, resembling a combination of the experimental surfaces proposed by Roscoe and Hvorslev, as shown in Roscoe et al. [78, 83].

The elastoplastic model with Drucker-Prager criterion and a cap offers improved predictive accuracy by considering the realistic compression behavior of soils, as discussed previously. Among the compression parameters, yielding stress emerges as the most important influential factor in determining the CBR value, as shown in Fig. 5 and Table 3. Figure 8 illustrates that significant portions of the material experience compression stress and remain in an elastic state, while substantial shear stresses are confined to small areas adjacent to the punch. Thus, in the FEM model, only specific elements experience stress paths that reach the shear criteria. These criteria depend on the friction angle and cohesion of the material. Therefore, it is possible to establish a correlation between the CBR value and Young's modulus and yielding stress. In addition, the CBR value is typically calculated at a strain of 2.5 mm. The present study proposes Eq. 16, which bears resemblance to the equation put forth by Mendoza & Caicedo [19] albeit with a nonlinear effect. This equation captures the material behavior based on the sensitivity analysis performed in Sect. "Sensitivity analysis of the input parameters predicting the CBR value".

$$CBR[\% ] = - 4 \times 10^{{ - 10}} p_{p} ^{{0.098}} E^{2} + \left( {0.002\;ln\;p_{p} - 0.00086} \right)E$$
(16)

Transportation geotechnics often involves correlations between the CBR value and the elastic modulus, as shown in Sect. "Literature review on factors and modeling that influence CBR value for coarse-grained soils". However, the paper also presents a variability analysis specifically focused on coarse-grained soils, examining the sensitivity of Young's modulus and the CBR value. The initial sensitivity analyses, depicted in Fig. 10, demonstrate that the CBR value increases with an increase in Young's modulus across all cases. Nonetheless, some studies have highlighted that the CBR value increment is not linear [6, 19, 84, 85]. This non-linear behavior is significant because many popular correlations, assume a linear relationship between Young's modulus and the CBR value [1,2,3]. Consequently, Fig. 11a provides comparisons between correlations derived from in pavement design guides, international standards, international recommendations, and elastoplastic simulations, revealing a significant amount of scattering in the CBR and Young's modulus correlations.

Fig. 10
figure 10

Capacity of elastoplastic models to represent CBR test

Full size image
Fig. 11
figure 11

a Comparison between practical correlations and CBR from simulations; b probability of underestimating values CBR from practical correlations

Full size image

In Fig. 11b, an academic exercise was done where the probability of obtaining higher CBR values (from the FEM simulations) was calculated for the analyzed correlations. This analysis assumed a normal probability distribution of the CBR values. Subsequently, the mean and standard deviation were obtained to calculate the probability of higher values for each correlation. From Fig. 11, it can be seen that the CBR values obtained from the correlations presented by Nielson [2], Heulelon & Klomp [1], MEPDG [4, 5], and CEBTP [3] have a probability of almost 100% of obtaining lesser values with respect to the FEM simulations for low CBR values (lower to 40%). The equation presented by Mendoza and Caicedo shows a less than 53% probability. The above changes with the increase in the expected value of the CBR with values of 68–82% by Nielson, 94.6% by the MEPDG, 98.8% by Heulelon & Klomp, and 99.7% by the CEBTP. In contrast, the correlation presented by Mendoza and Caicedo shows the opposite behavior by increasing the probability to 71% of underestimating the CBR value. So, the least conservative correlations to obtain the CBR value from Young's modulus are Mendoza and Caicedo, Nielson, MEPDG, Heulelon & Klomp, CEBTP. Still, in all cases, the probability of obtaining lower CBR values is less than 50%. It is clear that the analyzed correlations come from experimental work for different types of soils. As previously mentioned, all simulations conducted in this study were based on crushed hornfels rock [23]. The initial CBR results reported by Araya et al. [23] and subsequent calibration of parameters and variables allowed us to compute CBRs close to those obtained experimentally, as shown in Fig. 10.

Conclusions

The research used some of the more common constitutive models in geotechnical engineering to describe the mechanical behavior of coarse soil. These models use typical parameters used by geotechnics professionals with physical explanations. In addition, the simulations consider the intrinsic variability of the real soil. The simulations show that the increase or decrease in the CBR value is a function of the elasticity modulus, yielding stress and friction angle. Hence, an increase or decrease in CBR value can show a change in the previously mentioned parameters.

The simulations with parameter variability expand the knowledge of the shearing mechanisms, generated stresses, displacement fields, and load sharing when the CBR test is made. From these results, a physical explanation of test results can be drawn. FEM simulations showed the stress distribution zones in the soil during the CBR test. The more area is an elastic zone, followed by a compression zone and a high-shear stress zone. These zones can explain the importance of elastic modulus, yielding stress, and friction angle in the CBR value.

Monte Carlo simulations captured the variability effect of geotechnical parameters into CBR values for different coarse-grained soils. Hence, a basic linear correlation with the elastic modulus alone is insufficient. The CBR value also depends on the yielding stress; this correlation does not fit a linear relationship. Therefore, this paper introduced a novel approach to calculating the CBR value based on Young’s modulus and yielding stress, variables with greater weight according to the simulations. Regarding behavior with two variables, the CBR decreases when low-yield stresses were observed, exhibiting nonlinear behavior in response to changes in the elastic modulus. Conversely, high-yield stresses demonstrate nearly linear behavior as the elasticity modulus changes. In future research, the results can be validated with experimental tests influenced by factors such as particle fracturing, particle size, particle shape, and material composition of grains and materials mix, as shown in previous research. These factors can impact the elastic modulus and the yielding stress and collectively determine the CBR value.

A sensitivity analysis with the variability of geotechnical parameters was performed to compare the correlations proposed in the literature for practical design purposes. From these comparisons, determine the probability of underestimating values CBR was according to the correlation used to provide more accurate results. In most cases, these correlations can underestimate the behavior of the material.

Based on the results of this research, planned works involving more simulations of the new compounds to improve soil characteristics will be carried out in the near future. For example, biocatalysts, geopolymers, limes, recycled asphalt pavements, construction-Demolition Materials, and others. Based on the research presented by Lima et al. [27], Amaludin et al. [28], Putra et al. [29], Al-Obaydi et al. [30]. The results showed the important parameters that increase or decrease CBR value, and these parameters can be affected by the inclusion of this new material into the soils. As future results, equations that consider these new materials can be obtained. Simulators can also consider the effect of variability from the different proportions of the new material in the mix with the soil to optimize the design to improve the quality of materials. The impact can be using new technologies to improve the quality of materials in a more accessible form.