#### Abstract

The resilient modulus of subgrade is a design parameter of the pavement structure, which is significantly affected by the overlying load and traffic load. It is important to calculate the equivalent resilient modulus of the top surface of subgrade based on the nonuniform distribution of resilient modulus in subgrade. This paper takes the fully weathered granite soil as the research object. Firstly, the soil density of different layers of the subgrade structure is calculated by the degree of compaction of different subgrade layers. Secondly, the overlying load of each point in the subgrade is determined based on the quality of subgrade. Thirdly, the subprogram of the finite element software is compiled and redeveloped based on the elastic constitutive model, and the calculation method for the resilient modulus of each point in the subgrade under the traffic load is proposed when the convergence criterion is set up. Finally, according to the deflection equivalence of the elastic double layer and elastic half-space, the calculation and control methods for equivalent resilient modulus of the top surface of subgrade under nonuniform stress distribution are put forward, and the field verification tests are carried out. The results show that the error range between numerical calculation and field measurement of equivalent resilient modulus of subgrade is 10%. It shows that the calculation method for equivalent resilient modulus of subgrade proposed in this study is reasonable and effective. The equivalent resilient modulus of subgrade decreases with the increase of traffic load. And the resilient modulus of subgrade soil increases with the increase of subgrade depth. The resilient modulus of subgrade soil has a significant impact on the equivalent resilient modulus of subgrade after the overlaying gravel layer. The equivalent resilient modulus of the subgrade with the gravel layer increases with the increase of the resilient modulus of the subgrade soil.

#### 1. Introduction

The resilient modulus of all pavement layers, including subgrades, is one of the primary material properties as it is used in mechanistic pavement thickness design [1]. The resilient modulus (*M*_{r}) is widely applied and defined as the slope between the largest resilient strain and the largest deviatoric stress of dynamic triaxial tests [2]. Current studies suggested that the resilient behaviour of subgrade soil is affected by many conditions, such as the stress level and path, the moisture condition, the type of soil, and the soil structure [3, 4]. It is an acknowledgement that the stress state plays a key role in the *M*_{r} of subgrade soil [5]. For past decades, many researchers focused on the relationship between the *M*_{r} and stress states, and over a dozen of prediction models have been put forward with stress conditions [6]. Taking into account volume stress (*θ*), the *k*-*θ* model was originally mentioned by Seed et al. [2], but it is suitable only for the linear elastic behaviour. Witczak and Uzan introduced the octahedral shear stress (*τ*_{oct}) into the *k*-*θ* model to represent the shear effect between soils under the loading [7]. In 2004, the modified model was developed to overcome the indeterminate value problem of the Witczak model (*τ*_{oct} = 0) and applied into the structural design of highway subgrades by the Mechanistic-Empirical Pavement Design Guide (MEPDG) [8]. For multiple factors on subgrade soil, most of the predicting models were improved on the MEPDG model by embedding different stress terms and fitting coefficients [9].

Subgrade soil mainly withstands the impact of traffic load transferred from the road surface. As traffic load is transferred from top to bottom, stress decreases gradually [10, 11]. Thus, the impact of traffic load varies among soils of different depths. Alabi [12] and De Barros and Luco [13] calculated the vertical stress of subgrade under the action of vehicle load while studying the impact depth of vehicle load, showing that the vehicle load gradually decreases as the depth increases. Tang [14] and Han et al. [15] studied the mechanical behaviour of the dynamic load of subgrade soil and the response characteristic of dynamic stress while analyzing the law governing the impact on the dynamic stress of subgrade exerted by the value of dynamic vehicle load, frequency of dynamic load, different depths of subgrade, and number of dynamic load actions. It was also proved that the dynamic stress varies among different depths within the subgrade work area and gradually decreases as the depth increases [16].

The resilient modulus of subgrade is an important parameter to characterize the performance of subgrade [17]. However, many previous studies have also proved that the prediction model of resilient modulus of subgrade soil is highly dependent on stress states [18, 19]. The most famous prediction model of resilient modulus of subgrade is the three-parameter model proposed by Lytton, as shown in equation (1), which is ultimately included in the Mechanistic-Empirical Pavement Design Guide (MEPDG) [8]. Because of different stress states inside the subgrade, the *M*_{r} of subgrade soil is different at different points inside the subgrade structure. How to characterize the equivalent resilient modulus of subgrade is a very significant work in subgrade design.

The value of the *M*_{r} in the subgrade structure has a significant effect on the mechanical performance of the road surface. Because the resilient modulus of the soil at various positions of the subgrade structure is different under the dynamic loading, the equivalent resilient modulus of the top surface of subgrade is usually used as the design parameter of subgrade and pavement structure. In the last decades, many researchers have concluded that the decrease of the equivalent *M*_{r} of subgrade can lead to an increase in surface deflection [20]. Sadrossadat et al. [21] and Gu et al. [22], who studied the impact of the *M*_{r} of subgrade on the structural performance and mechanical response of the road surface, argued that the *M*_{r} of subgrade exerts an important influence on the service life of the pavement structure. The change of the climatic environment in four seasons is accompanied by a change in the water content of subgrade soil. Thus, the stiffness of the subgrade decreases significantly and the deformation increases continuously because of the effect of the climatic environment and cumulative load [23]. Some researchers observed that, on placing the unbonded granular layer on the top of the embankment, the *M*_{r} of subgrade construction will increase [24, 25]. However, the traffic load associated with the elastic layered system is always regarded as a static load, failing to reasonably reveal the regulation and control of the equivalent *M*_{r} of subgrade. The impact of dynamic traffic load cannot be dismissed because of the nonlinearity and anisotropy of the unbonded granular layer and nonuniform spatial distribution of the *M*_{r} inside the subgrade [26, 27].

To sum up, the resilient modulus of subgrade soil strongly depends on stress; however, the spatial distribution of the resilient modulus inside the subgrade is nonuniform because of the nonuniform distribution of static load and dynamic stress inside the subgrade; few researchers focus their research on this topic, which leads to the lack of a calculation method for the equivalent resilient modulus of subgrade based on nonuniform spatial distribution of stress, let alone a regulation and control method in this regard. Therefore, taking into account the nonuniform stress field in the subgrade, it is essential for the pavement design to calculate accurately the equivalent resilient modulus of the subgrade. Meanwhile, it is reasonable that increasing the stiffness of the roadbed can effectively control the deformation and fatigue life of the pavement structure.

#### 2. Calculation Method for Equivalent Resilient Modulus of Subgrade Based on Nonuniform Distribution of Stress

The dynamic load and overlying load on the subgrade soil vary among different depths within the subgrade work area. Therefore, the state of confining pressure (overlying load) and the state of dynamic load (traffic load) vary among different points inside the subgrade. These two stress states exert impacts on the resilient modulus of the subgrade soil, which is also the main reason for the difference of resilient modulus among various points inside the subgrade soil and results in the nonuniform distribution of the resilient modulus inside the subgrade soil. Before the calculation of the equivalent resilient modulus of subgrade, it is necessary to first accurately calculate the resilient modulus of various points inside the subgrade under the nonuniform distribution state and then calculate the equivalent resilient modulus of the top surface of subgrade.

##### 2.1. Determination of Soil Density at Different Layers of Subgrade

Before the determination of the overlying load at various points of the subgrade, it is necessary to calculate the soil density at different layers of the subgrade. Generally speaking, subgrade soil is filled at certain compactness while the optimum water content of the soil is reached [28]. It includes three layers, which are the roadbed, upper embankment, and lower embankment. The compactness of subgrade varies at different positions of the subgrade. By the current highway subgrade design specification of China [29], the wet density of soil at various positions of subgrade is shown in equations (2)–(4).

Roadbed:

Upper embankment:

Lower embankment:where is the maximum dry density of the soil, *θ*_{iw} is the volumetric water content of each point in the subgrade, and *n* is the number of nodes of each layer of the subgrade in the numerical calculation model.

##### 2.2. Determination of Overlying Load at Various Points inside Subgrade

To determine the overlying load on various nodes inside the subgrade is a relatively complex process. Formulas for calculating the overlying load on various nodes of the roadbed, upper embankment, and lower embankment are as follows.

Roadbed:

Upper embankment:

Lower embankment:where *y*_{i} is the vertical coordinate value of any point *i* inside the subgrade (unit: m); is the maximum value of the vertical coordinate in all node coordinates of the subgrade model; is the gravitational acceleration, taking 10 m/s^{2}; *ρ*_{0} is the density of the pavement structural layer, taking 2500 kg/m^{2}; *ρ*_{1} is the density of the roadbed; *ρ*_{2} is the density of the upper embankment; *ρ*_{3} is the density of the lower embankment; and *h* is the thickness of the pavement structural layer.

##### 2.3. Determination of Resilient Modulus of Various Points inside Subgrade

###### 2.3.1. Development of User Subprogram in ABAQUS

A subgrade structure model is created through the ABAQUS software. Subgrade soil is subject to the isotropic linearity-elasticity constitutive relation, with the stiffness matrix shown as

As the subgrade withstands traffic load transferred to the road surface, the dynamic load varies among soil of different depths, resulting in the nonuniform distribution of the modulus field inside the subgrade. For finite element calculation, each grid cell can be regarded as an elastomer; that is to say, the resilient modulus of each grid cell is a fixed value. To determine the resilient modulus of each grid cell, it is necessary to use the UMAT program for secondary development in ABAQUS. The specific method is given below.

The initial modulus (taken as 60 MPa herein) is set for each grid cell in the subgrade structure, the standard axle load is given, and then the first trial calculation is undertaken, whereby a dynamic stress *q*_{icyc} is assigned to each node. Based on the authors’ previous study [6], the relation between octahedral shear stress and deviatoric stress is . Therefore, the octahedral shear stress of each node in the subgrade structure can be obtained, i.e., . And the *M*_{r} of each node (grid cell) can be obtained, as follows:

Many researches have studied the nonlinear convergence of the elastic modulus. The tangent stiffness method was applied to solve the nonlinear stress-strain behaviour by updating the tangent stiffness at each iteration until the stress increment converges [1]. Young’s modulus *E* (Young’s modulus is equal to the resilient modulus in a pure elastomer) in equation (8) is replaced by *M*_{ir} corresponding to each node to obtain a new modulus of each node. The new modulus of each node is used to update the previous initial value of modulus, and the next dynamic load is then applied. This process is repeated. In Figure 1, the process of the UMAT subroutine was simplified to ensure that the style of the flow chart is consistent. In fact, the update of the *M*_{r} not only occurs in a single node but also occurs in all the FEM system.

The modulus of each node inside the subgrade varies under different dynamic stresses. When the error in the modulus between two calculations for a definite node is within 0.5%, and the accumulated error in the modulus between two calculations for all the nodes of the subgrade structure is within 5%, the accuracy of the modulus of each grid cell inside the subgrade is considered to meet the project requirements. The finite element calculation is concluded at this point. The convergence criteria for the modulus field inside the subgrade are as follows:where *M*_{r}^{i} is the output modulus of the cycle *i*; *M*_{r}^{i−1} is the output modulus of the cycle *i*−1; error_{i} is the error of the *M*_{r} of each node (grid); error_{c} is the cumulative error of the entire model; and *n* is the number of cycles.

The flow chart for the calculation of the nonuniform modulus field inside the subgrade is shown in Figure 1.

###### 2.3.2. Deflection at the Top of Subgrade under Uniform Modulus

For the calculation of the deflection of the top surface of subgrade under the uniform distribution of resilient modulus of subgrade, the above-mentioned finite element calculation model and its relevant parameters are adopted. The resilient modulus field of subgrade does not vary with the stress state in this model. To numerically calculate the deflection (deformation) of the top surface of the subgrade, the dynamic loading applied on the top surface of the subgrade is half-sinusoid cyclic loading with a stress amplitude of 0.7 MPa, a loading cycle of 0.2 s, an interval time of 0.8 s, and a loading area of 30 cm diameter circle. The cycle period is one second. The loading conditions mentioned above are used to calculate the resilient modulus of subgrade soil at 20∼120 MPa, and the deflection at the central point of stress applied on the top surface of the subgrade can be obtained. The relation between the resilient modulus of subgrade soil and the deflection is shown in Figure 2.

**(a)**

**(b)**

#### 3. Site Verification of Equivalent Resilient Modulus of Subgrade Based on Nonuniform Distribution of Stress

The site verification with a portable falling weight deflectometer (PFWD) is carried out on the top surface of the upper embankment (subsoil) of K116 + 980∼K117 + 180 (left part) in Lot C of Guangzhou–Foshan–Zhaoqing Highway. The PFWD site tests are also conducted on the top surface of unscreened gravel layers with different thicknesses laid on the upper embankment (subsoil).

##### 3.1. Testing Program

In this study, the PFWD site test section is 200 m in length. Firstly, two testing lines are placed on the top surface of the upper embankment along the vehicle wheelmark of the carriageway, and the horizontal distance between the two test lines is 2.0 m. The PFWD site tests are conducted on each testing line at a spacing of 20 m. The deflection of each test point on the top surface of the subgrade can be obtained by the PFWD. There are 24 test points in total. Secondly, after each PFWD test on the top surface of the upper embankment (subsoil) was carried out, the unscreened gravel layer of 10 cm thickness is laid on the upper embankment. And the PFWD site tests are also done in the same positions on the surface of the unscreened gravel layer of 10 cm thickness. Thirdly, the above experimental methods had been reused on the top surface of the unscreened gravel layer of 20 cm thickness and 30 cm thickness. There are 72 test points in total, and the stress-time curves of the dynamic loading can be obtained. And also the measured deflection of each test point in each layer can be gained.

##### 3.2. Material Parameters for Subgrade Soil

The filling subgrade soil is granite eluvial soil (sandy clay with a low liquid limit). After the PFWD tests on the top surface of the subgrade (subsoil) are implemented, the unscreened gravel is placed on the top surface of the subgrade and compacted to the layer with a thickness of 10 cm, and the compactness of the unscreened gravel layer is not less than 96%. After that, the test points are located by the electronic total station. The site tests of the PFWD are conducted in the same positions when the unscreened gravel layers are laid on the top surface of subgrade of 10 cm thickness, 20 cm thickness, and 30 cm thickness. More than 5 PFWD tests are carried out at each test point, and the average value of tests is taken as the measured value for every test point.

##### 3.3. Finite Element Modelling of Subgrade Structure of Test Section

Finite element model of subgrade is undertaken according to the size and parameters of the test section in Figure 2(a). In the ABAQUS model, all mesh and nodes are divided by a series of 8-node quadratic reduced integration elements. The element mesh is inherited from the attribute of the humidity calculation model in the GeoStudio/Vadose software, where a node is divided at 0.5 m intervals. The subsidence of the simulated surface of the road was caused by the vehicle’s axial load, without considering the influence of self-gravity. The bottom and the side of the model were applied to the constraints.

To study the impact of unscreened gravel on the equivalent resilient modulus of subgrade for different thicknesses of unscreened gravel layers, finite element models with different thicknesses of unscreened gravel layers are created. The method mentioned above for the finite element calculation of equivalent resilient modulus of subgrade is used for numerical calculation of the standard half-sinusoid loading curve with a stress amplitude of 100 kPa and a dynamic loading cycle of 0.2 s.

##### 3.4. Comparison and Verification of Calculated and Measured Values for Top Surfaces of Subgrade and Unscreened Gravel Layers

To investigate the difference between the calculated values and the PFWD-measured values of the deflection of the top surface of subgrade, a half-sinusoid loading curve under standard axle load is obtained according to the above-mentioned method. The calculated deflection value of the top surface of the subgrade is 152.4 (0.01 mm), measured by the PFWD tests. The test results are shown in Tables 1 and 2.

After the deflection of the top surface of subgrade soil in the test section is measured, unscreened gravel is placed on the top surface layer by layer; there are three layers in total, and each layer is 10 cm thick. The PFWD test is conducted on each layer. The coordinates of the test point correspond to those of the test point on the top surface of the subsoil.

When each unscreened gravel layer (10 cm, 20 cm, and 30 cm in thickness) is placed on the top surface of the subgrade, the deflection of the top surface of the subgrade is calculated. The material resilient modulus of the unscreened gravel layers is 200 MPa, which is obtained through the laboratory test [30].

It can be clearly found after comprehensive analysis of Tables 1 and 2 that, before unscreened gravel layers are placed, the average deviation value between the calculated deflection and the PFWD-measured deflection at the central position on the carriageway of the top surface of the subgrade is 12.2 (0.01 mm), with an error rate of 7.41%. Moreover, this difference at the position on the wheelmark belt of the carriageway is 11.3 (0.01 mm), with an error rate of 6.93%; the average deviation value of structural modulus is 4.60 MPa, with an error rate of 7.17%. Then, through conversion with the data in Figure 2(b), the average difference of structural modulus is 5.27 MPa, with an error rate of 8.31%. After 10 cm, 20 cm, and 30 cm unscreened gravel layers are placed, the average deviation value between the calculated deflection and the PFWD-measured deflection at the central position on the carriageway of the top surface of the subgrade is 11.3, 10.5, and 10.0 (0.01 mm), respectively, with an error rate of 7.83%, 7.92%, and 8.09%, respectively.

From Tables 1 and 2, it can be seen that, before unscreened gravel layers are placed, the average deviation value between the calculated deflection and the PFWD-measured deflection at the position on the wheelmark belt of the carriageway is 11.3 (0.01 mm), with an error rate of 6.93%; the average deviation value of structural modulus is 4.60 MPa, with an error rate of 7.17%.

#### 4. Control Method for Equivalent Resilient Modulus of Subgrade

Considering that many previous researchers paid more attention to the resilient modulus of the overlay or calculated the equivalent resilient modulus of the gravel layers under static loading, and the research results on the equivalent modulus of the gravel layer top under dynamic loading are few, especially considering the nonuniform distribution of the resilient modulus of subgrade soil, there are few related in-depth studies. One of the purposes of this paper is to provide a quantitative variation rule of equivalent resilient modulus for road engineering researchers under the condition of nonuniform distribution of resilient modulus in subgrade.

##### 4.1. Impact of Different Subsoil Modulus Values and Layer Thicknesses on Structural Modulus

To study the impact of different material modulus values and layer thicknesses on the equivalent resilient modulus of subgrade under different subsoil modulus conditions, a finite element calculation model is created for different material modulus values and layer thicknesses, based on the cross-sectional size and subsoil material parameters of the test section of Nanchang–Zhangshu Highway. The working conditions for calculation are shown in Table 3. The equivalent resilient modulus of the surface of the subgrade structure is calculated. The results of calculation are shown in Figure 3.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

According to Figure 3, before unscreened gravel layers are placed on the subsoil, the subsoil modulus has a significant impact on its overall equivalent resilient modulus. At the layer thickness of 30 cm and the layer modulus of 300 MPa, as the subsoil modulus increases from 20 MPa to 120 MPa, its equivalent resilient modulus increases from 32 MPa to 143 MPa, about 4.5 times. This shows that the subsoil modulus is an important factor in affecting its equivalent resilient modulus, and such layers can effectively increase the equivalent resilient modulus of the subgrade. With the fixed subsoil modulus, the equivalent resilient modulus of the subgrade increases as the modulus of such layers increases. When the *M*_{r} of the subsoil is 60 MPa and layer thickness is 30 cm, as the layer modulus increases from 100 MPa to 500 MPa, the equivalent resilient modulus of the subgrade also increases from 62 MPa to 82 MPa (up about 30%). With the fixed subsoil modulus and fixed layer modulus (larger than subsoil modulus), the equivalent resilient modulus of the subgrade increases as the layer thickness increases. At the subsoil modulus of 60 MPa and the layer modulus of 300 MPa, as the layer thickness increases from 0 cm to 30 cm, the equivalent resilient modulus of the subgrade increases from 59 MPa to 66 MPa, up about 15%. According to Figure 3 and the above-mentioned analysis, when the subsoil modulus is relatively low, layers with a larger material modulus can significantly increase the overall equivalent resilient modulus of the subgrade; when the subsoil modulus is relatively high, layers with a higher material modulus is less effective in increasing the overall structural modulus of the subgrade.

##### 4.2. Law Governing Change in Structural Modulus of Subgrade under Different Material Modulus Values and Layer Thicknesses

To study the influence of different material resilient modulus values and layer thicknesses on the equivalent resilient modulus of subgrade under different subgrade soil modulus conditions, the same calculation model and parameters in Table 3 are used. For the finite element numerical calculations under different soil modulus values and thicknesses, the calculation conditions are shown in Table 3 and the results of the equivalent resilient modulus of the subgrade are shown in Figure 4.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

According to Figure 4, when unscreened gravel layers are placed on the subsoil, the material modulus of the layers exerts significant influences on the overall equivalent resilient modulus of the subgrade. At the subsoil modulus of 40 MPa and the layer material modulus of 300 MPa, as the layer thickness increases from 0 cm to 30 cm, the equivalent resilient modulus of the subgrade increases from 40 MPa to 50 MPa, up about 25%.

In conclusion, layers with a certain thickness and higher material modulus on the top surface of subgrade (top surface of subsoil) can effectively increase the overall equivalent resilient modulus of the subgrade. Though layers cannot increase the equivalent resilient modulus of the subsoil, we can design corresponding layers according to the attenuation of the subsoil modulus during road operation, propose a regulation and control method for the equivalent resilient modulus of the subgrade, increase the overall equivalent resilient modulus of the subgrade, and reduce the attenuation rate of the structural modulus of the subgrade during its service life so as to ensure the stability of the structural modulus of the subgrade throughout the road operation period.

#### 5. Conclusions

(1)Based on the comparison between the PFWD-measured deflection and the calculated deflection of the top surface of the subgrade, the error values between calculated and measured results of the equivalent resilient modulus of the subgrade are within 10%. This demonstrates that the calculation method for the structural modulus of subgrade proposed in this study is reasonable and effective.(2)The equivalent resilient modulus of the subgrade reduces as the traffic load increases; the material modulus of subgrade soil inside the subgrade increases as the subgrade depth increases. When the load is less than 2.0 times the standard load, the equivalent resilient modulus of the subgrade reduces more significantly as the load increases; within the vertical depth of 1.5 m, the resilient modulus of the material of the subgrade soil increases more significantly as the depth increases, at about 60%.(3)The resilient modulus of subgrade soil has important impacts on the overall equivalent resilient modulus of the subgrade provided by laid layers. In response to the fact that the resilient modulus of subgrade gradually reduces with climatic changes, we have proposed a control method for the equivalent resilient modulus of the subgrade; that is, the gravel layers with a certain thickness, higher modulus, and good water stability are placed on the top surface of the subgrade to increase the equivalent resilient modulus of the top surface of the subgrade. And it can ensure the stability of the resilient modulus of subgrade throughout the road operation period.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Disclosure

The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### Conflicts of Interest

The authors declare no conflicts of interest.

#### Acknowledgments

The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China (51878078, 51608057, and 51908562), Hunan Provincial Innovation Foundation for Postgraduate (CX2018B521), Open Research Fund of Science and Technology Innovation Platform of State Engineering Laboratory of Highway Maintenance Technology, and Changsha University of Science & Technology (kfj150103). This study was also supported by the State Bureau of Forestry 948 Project (2015-4-38), the Youth Scientific Research Foundation, Central South University of Forestry and Technology (QJ2018008B), and the Open Fund of Engineering Research Center of Catastrophic Prophylaxis and Treatment of Road & Traffic Safety of Ministry of Education (Changsha University of Science & Technology) (kfj170405).