# Back-Calculation of Soil Modulus from PFWD Based on a Viscoelastic Model

Authors: Rui Zhang; Tianzeng Ren; Md Ashrafuzzaman Khan; Yusheng Teng; Jianlong Zheng
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

To improve the accuracy of back-calculation of soil modulus using the portable falling weight deflectometer (PFWD), a viscoelastic method (VEM) overcoming the limitations of the conventional linear elastic method (LEM) was proposed. A quasi-static dynamic analysis technique of Laplace transformation and a modified Gauss–Newton optimization algorithm were adopted in the proposed method. The back-calculation program was compiled with MATLAB. The effectiveness of the proposed method was verified with the in situ plate load test (PLT) conducted on a highway embankment. In situ test results showed that a time lag existed between the peaks of deflection and load, and load-deflection curves were nonlinear, which indicated the viscoelastic nature of the soil. The back-calculated modulus using the LEM and VEM was higher than that using the PLT. In the case of low stiffness soil, the average error of back-calculation using the LEM and VEM was 53.1% and 14.8%, respectively. However, for stiffer soil, the average error of back-calculation using the LEM and VEM method was 12.4% and 4.3%, respectively. Moduli of back-calculation using LEM and VEM methods were used to perform flexible pavement analysis, which showed that with an 8% reduction of modulus, the pavement service life reduced by 25%. More accurate estimation of modulus can save maintenance cost in the future.


---

## Body

## 1. Introduction

The Mechanistic-Empirical Pavement Design Guide and other existing pavement design guides use elastic modulus (E) as the primary input parameter for subgrades. The soil modulus is typically determined by different laboratory tests or in situ nondestructive tests. The falling weight deflectometer (FWD) and plate load test (PLT) can be used to determine the stiffness of pavement, but these devices have some limitations when they are used during the construction of pavement subgrade. The FWD is not always applicable to the project site. And the PLT takes several hours so that it can slow down the construction process. These limitations can be overcome by conducting a test with a portable falling weight deflectometer (PFWD), a nondestructive and cost-effective tool used to determinate the soil modulus [1–6]. The soil modulus can be determined by the deflection profile and the back-calculation method.Over the past decade, the PFWD as a portable and cost-effective tool for determining the stiffness of foundation soil has gained much popularity [1–3, 7, 8]. Several static and dynamic back-calculation programs use the deflection profile to perform back-calculation analysis for determining the in situ elastic modulus. However, the modulus obtained from these programs may vary because of different assumptions, interactive technique calculation schemes, and material models [1, 4, 7, 9–14]. Stress and strain calculated from the available methods differ from the actual values in the pavement material and depend on the material model adopted for the pavement analysis [12, 15–17]. Although the base/subgrade is considered as a continuous, homogeneous, linear elastic materials in most of studies; in reality, they exhibit noncontinuous, nonhomogenous, anisotropic, viscous, plastic, and viscoelastic behaviors.Based on the multilayered elastic theory, ASTM D5858-96 [18] provides a guideline for back-calculating the elastic modulus by deflections. Most researchers [1, 9, 12, 17, 19] used the linear elastic material model in their studies, as no additional guidance is available from the American Association of State Highway and Transportation Officials or the American Society of Testing Materials.Modulus can be back-calculated by the available computer programs (e.g., MODULUS, MODCOMP, WESDEF, WESNET, MICHBACK 1.0, FWD-DYN, ELMOD, EVERCALC, WESDEF, CHEVDEF, ELSDEF, FDEDD1, SEARCH, FPEDDI, ISSEM4, PDADAL, and BISDEF). Most of these programs are based on the linear elastic layer theory, but their accuracy depends on the goodness of fit between the computed and measured deflections [15]. Assuming that the surface load is static and uniformly distributed, stress and strain are computed by these programs. Layers are linear elastic, isotropic, and homogeneous.The PFWD test is dynamic. It has been observed that the linear elastic material model cannot accurately predict the response under dynamic load. Asli et al. [20] used a viscoelastic soil model to calculate the soil modulus. The plate/soil system is represented by a single degree of freedom (SDOF) system. A least square algorithm was used as an error minimization technique. Elastic stiffness, the ratio of applied static load to measured deflection, can be estimated by the PLT. As the PFWD test is dynamic, the dynamic response of ground is affected by inertia and damping properties of media. Asli et al. [20] adopted a dynamic equilibrium equation in which the magnitude of applied load depends on acceleration, velocity, deformation, damping coefficient and stiffness of elastic medium, and applied load on the top of it. Based on the least-squares method, the elastic stiffness, damping coefficient, and equivalent mass were determined by the minimization technique, but the time parameter was omitted. Because of the lack of a time parameter, based on this method, the actual viscoelastic behavior could not be predicted. It was also assumed that the subgrade was considered as the Winkler foundation, a system of several individual springs. There are no internal connections between the springs, and the deformation of a certain point only depends on the pressure acting on the point. However, there is a direct transverse relationship between the soil particles, which restricts each other in the transverse direction. When the vertical load is applied to a certain point, it will deform vertically along with the adjacent soil particles and make a deflection basin instead of a single point deflection. Therefore, the soil medium considered as the Winkler spring system cannot predict the actual deflection shape.There are two methods to analyze the dynamic response of viscoelastic system: quasi-static and dynamic methods. Amir and Javad [21] developed a dynamic analysis framework, which was used for any explicit/implicit time integral. Hudsadin and Yasothorn [22] used a dynamic half-space model for dynamic back-calculation of engineering soil properties based on the FWD data. In dynamic back-calculation, the structure of equation is complex, which is simplified by a numerical algorithm, that is, a time-consuming calculation method. Ameri et al. [23] adopted a layer-wise ﬁnite element equation for quasi-static analysis of laminated structures with embedded viscoelastic material whose constitutive behavior is represented by the Prony series. Ahmet et al. [24] studied the quasi-static behavior of viscoelastic Kirchhoff plates by the finite element method in transformed Laplace–Carsonspace. Compared with the dynamic method, the quasi-static method is simpler and less time-consuming.Based on a viscoelastic model, a back-calculation method which incorporated a quasi-static dynamic response analysis method with elastic half-space theory was proposed first. The proposed method was then verified by field test results. Finally, a case study was conducted to analyze the importance of choosing an accurate soil modulus in pavement design.


## 2. Field Studies

The Hainan Wanning-Yangpu Highway is 163 km long. Approximately 2.4 million m3 high plasticity soils were excavated and used as filling materials for the construction of embankment in many sections of highway. The section between K110 + 180 and K110 + 280 was selected as a trial section in this study.As shown in Figure1(a), the cross section of embankment was divided into two parts in the vertical direction. One is 1.5 m high on the top of upper embankment, and the other is at the lower embankment. High plasticity soil was used to construct the lower embankment. The basic properties of soil in the trial section were determined by the laboratory tests. The natural moisture content, contents of fine-grained particle, liquid limit, and plastic limit were 29.1%, 60.3%, 54.0%, and 28.5%, respectively. The optimum moisture content and the maximum dry density determined by the standard compaction method (ASTM D698-12e2 [25]) were 17.3% and 1.76 g/cm3, respectively. The high plasticity soil layers were compacted at a moisture content of 19% and its maximum dry density of 95%. The slopes of both sides were covered with low plasticity soil to prevent surficial slope failure. The upper embankment was constructed with high plasticity soil mixed with 4% cement. The optimum moisture content and the maximum dry density determined by the modified compaction method (ASTM D1557 [26]) were 15.6% and 1.82 g·cm−3, respectively. The cement-treated soil layers were compacted at a moisture content of 17.3% and its maximum dry density of 93%.Figure 1
Embankment structure diagram and layout of measuring points. (a) Cross section of embankment. (b) Top of the lower embankment (elevation is 5.0 m from the ground level). (c) Top of the upper embankment (elevation is 6.5 m from the ground level).
(a)
(b)
(c)The modulus of high plasticity soil and cement-treated soil in compaction state was determined by the dynamic triaxial test (AASHTO T307-99 [27]) (Figure 2). Three samples were selected in this test. The average moduli of high plasticity soil and cement-treated soil were 19.6 MPa and 201.2 MPa, respectively. The coefficient of variation is about 5%.Figure 2
Dynamic triaxial test.After the construction of lower embankment, the PFWD and PLT (Figure3) were conducted in 10 specific locations (L1 to L10) and the spacing between two locations nearby was kept to 10 m, as shown in Figure 1(b). Similarly, after the completion of upper embankment, the field tests were carried out on 10 locations (U1 to U10), as shown in Figure 1(c).Figure 3
Field tests on the trial embankment of Wanning-Yangpu highway. (a) PFWD test. (b) PLT.
(a)
(b)The primary parameters of PFWD equipment are as follows: the load plate radius is 0.15 m; the load range is 1–15 kN; the load pulse shape is half-sine wave; the load sensor relative error is less than 1%; the absolute error is ±0.1 kN; the effective test accuracy of deflection sensor is 1μm; the maximum test range is 2200 μm; and the data acquisition period is 0∼60 ms. The PFWD test methods and procedures were carried out according to ASTM E2583-07 [28]. The PLT was conducted in accordance with relevant provisions of the ASTM D1195 [29].The load-deflection and time-history curves obtained from these locations were similar in nature, and the data obtained from the L4 and U4 points were selected for further analysis. The field test results obtained from the top of lower and upper embankments are shown in Figures4 and 5, respectively.Figure 4
PFWD results for L4. (a) Load and deflection time-history curve. (b) Load-deflection curve.
(a)
(b)Figure 5
PFWD results for U4. (a) Load and deflection time-history curve. (b) Load-defection curve.
(a)
(b)The time lag between the peaks of deflection and load for the low stiffness soil can be observed from Figure4(a). Due to the viscoelastic behavior of high plasticity soil, the peak deformation was 3 ms later than the peak of load. When the applied load reaches the peak, the corresponding deformation is lower than the peak value; When the deformation reaches the peak, the corresponding load is lower than the maximum load. The magnitude of load curve started to decrease after reaching the peak, but the magnitude of deformation curve still increased to the peak during the lag time. Nonlinear behavior and hysteresis loops were observed in the load-deflection curve, as shown in Figure 4(b). Therefore, the deformation hysteresis effect of high plasticity soil under the impact load is significant.The PFWD and PLT were carried out on the upper embankment surface of K110 + 180∼K110 + 280 section to study the applicability of VEM to the stiffer soil embankment test. The PFWD test results at U4 are shown in Figures5(a) and 5(b). Compared with Figures 4(a) and 4(b), due to the increase of soil modulus, the deformation hysteresis effect and load-deflection hysteresis loop were significantly reduced, but they still existed, thus causing an error in the LEM back-calculation.In the existing conventional method, it is assumed that to back-calculate the elastic modulus by the PFWD, the relationship between load and deformation is approximately linear. The Boussinesq solution of elastic half-space in equation (1) was used to calculate the soil modulus:(1)E=πR2p1−μ2w,where E is the soil modulus value (MPa), p is the peak load plate pressure (MPa), R is the radius of the load plate (m), Poisson’s ratio μ is 0.35, and w is the peak vertical deﬂection (mm).For the top of lower embankment (L1–L10), the comparison of the LEM modulus with the measured modulus are shown in Figure6(a). The maximum, minimum, and average values of the LEM and measured moduli are 27.5 MPa and 18.4 MPa, 24.5 MPa and 15.1 MPa, 26.0 MPa and 17.0 MPa, respectively. For the lower embankment, the average absolute and relative errors are 9.0 MPa and 53.1%, respectively. The measured and back-calculation results of all 10 measuring points on the top of upper embankment (U1–U10) are shown in Figure 6(b). The maximum, minimum, and average values of the LEM and measured moduli are 238.1 MPa and 213.2 MPa, 208.2 MPa and 184.5 MPa, 221.8 MPa and 197.3 MPa, respectively. For the upper embankment, the average absolute error is 24.5 MPa, and the average relative error is 12.4%.Figure 6
Comparison of back-calculated moduli using LEM and measured moduli using PLT. (a) Results of L1–L10. (b) Results of U1–U10.
(a)
(b)The loading frequencies of the PFWD, dynamic triaxial and PLT are about 4 Hz, 1 Hz, and 0 Hz, respectively. When the soil is loaded at low frequency, the pulse duration is longer and the behavior of soil is mainly controlled by its quasi-static stiffness. However, when the pulse duration is short and the frequency increases, the soil behavior is mainly controlled by inertia, which weakens the dynamic response amplitude of soil. Only the peak values of load and vertical deflection were considered in equation (1), which generated a higher value of elastic modulus than the modulus obtained from the dynamic triaxial and PLT.


## 3. Proposed Moduli Back-Calculation Method Based on a Viscoelastic Model

The LEM is not accurate. Therefore, it is necessary to establish a simple method for back-calculation modulus that can better reflect the nonlinear dynamic response characteristic of high plasticity soil embankments under an impact load.


### 3.1. Selected Viscoelastic Theoretical Model

It is well known that under a sinusoidal dynamic load, the strain of viscoelastic material changes periodically, and the strain response lags behind the stress response. The viscoelastic theoretical model has been commonly applied to dynamic response analysis of asphalt pavement under the traffic load. The deformation hysteresis effect of a high plasticity soil embankment under impact load is similar to the mechanical behavior of viscoelastic materials mentioned earlier. Therefore, based on the viscoelastic constitutive model, a high plasticity soil modulus back-calculation method was established in this research.The Kelvin and Maxwell viscoelastic models are relatively simple, whereas the three-parameter solid model, Burgers model, and generalized Maxwell model are more complex. The Maxwell model is a series of spring-dampers that are primarily used to simulate the stress relaxation process of viscoelastic materials. The three-parameter solid model, Burgers model, and a generalized Maxwell model are constructed by a spring-damper parallel and series. The calculated parameters were incorporated to better describe the viscoelastic and residual deformation properties of material under dynamic load. However, the multiparameter model greatly increases the difficulty and iteration time of the back-calculation, but the results may not converge. In the PFWD test, the loading and unloading time of drop hammer impact is very short so that the plastic deformation cannot occur in time. The vertical deformation measured by PFWD is mainly rebound deformation. However, only the recoverable deformation was used to determine the elastic modulus, and the permanent deformation part was eliminated by the repeated dropping of the hammer. Therefore, it is not necessary to consider the residual deformation or adopt the complicated multiparameter viscoelastic model to calculate the high plasticity soil modulus.In this study, the Kelvin model was selected, and its concise structure was conducive to an efficient iteration of subsequent MATLAB program. The mechanical response of subgrade can be described by this model to meet the accuracy requirements.The Kelvin model, a viscoelastic model, was composed of a spring-damper in parallel. It is simple in form and easy to iterate efficiently. Its mathematical description is as follows:(2)σ=Eε+ηdεdt,where σ is the stress, ε is the strain, E is the elastic modulus, and η is the viscosity coefficient.The load-time curve obtained from the PFWD impact load can be considered as a semisine curve. Assuming that the mathematical expression of curve isσ¯sinωt, equation (3) was obtained by substituting σ¯sinωt into equation (2) and ignoring the rapid attenuation index items associated with time t:(3)σσ¯2−2cosδσσ¯εε¯+εε¯2=sin2δ,where σ¯ is the stress amplitude, ε¯ is the strain amplitude, δ is the phase difference, and tanδ=ηω/E.According to equation (3), the stress-strain curve under the action of loading and unloading conditions can be obtained, as shown in Figure 7. It could be seen that when the viscosity coefficient was introduced, the loading stress first reached the maximum value, and the corresponding strain did not reach the maximum. However, the maximum strain was observed during the unloading process due to the hysteresis. The stress-strain curve is very similar to the PFWD-measured load-deflection curve, as shown in Figure 4(b). Therefore, the Kelvin model was adopted to analyze the dynamic response of high plasticity soil under impact load.Figure 7
Stress-strain curve of Kelvin model.


### 3.2. Quasi-Static Back-Calculation Method

The dynamic and quasi-static methods were used for simulating dynamic response analysis of elastic systems under dynamic loading.In the dynamic back-calculation [20], combined the motion differential equation of axisymmetric space problem in the cylindrical coordinate system with the physical equation of stress expressed by displacement, the Navier dynamic equation can be obtained. The viscoelastic theory model was introduced into the Navier dynamic equation, and the explicit solution of vertical deflection was derived by Laplace and Hankel back-calculation. The complex mixed-boundary-value problem was also considered. Since the structure of implicit partial differential equations was complicated in back-calculation, time-consuming numerical algorithms were adopted to simplify it.For this research, the quasi-static back-calculation method is to take the explicit solution of vertical deflection on the surface of the elastic half-space body under the static load as Laplace transformation and to replace the elastic parameters in the equation with viscoelastic operators. In this method, the dynamic load data were used to calculate the static parameters. Subsequent back-calculation results verified that the quasi-static method met the requirements for back-calculation of high plasticity soil embankments. Therefore, the modulus calculated by this method and the PFWD dynamic load data described in this study is static modulus, which is consistent with the static modulus measured by the PLT.The vertical deflection on the surface of elastic half-space body under the action of an axial symmetric load was calculated as follows:(4)wr=21−μ2E∫0∞J0ξrp¯ξdξ,where J0(x) is the zero-order Bessel function, p¯ξ is the zero-order Hankel integral transformation of the load for PFWD rigid bearing plate, p¯ξ=pRsinξR/2ξ; p is the load plate pressure (kN⋅m−2), R is the radius of the load plate (m), and r is the radial distance from the load center of calculation point (m).After the Laplace transformation of equation (4) was performed, the viscoelastic operator was used to replace the elastic modulus, and the PFWD load function p¯ξ was substituted into the equation for inverse Laplace transformation. Since the Poisson ratio, which changes less, has little influence on the calculation result, it was not considered in the calculations. The vertical explicit deflection equation of viscoelastic half-space surface under the action of the PFWD rigid loading plate was obtained as follows:(5)wr,t=L−1w¯r,s=1−μ2Fπ⋅Es⋅R1−e−E/ηt∫0∞J0ξrsinξRξdξ,where E(s) = E + ηs is the viscoelastic operator; F is the concentrated force generated by the PFWD (KN) acting on the surface of the measuring point, which is obtained by the load time-history curve; R and r have the same meaning as equation (4); and s is the Laplace transformation parameter of time t.In equation (5), it is shown that when t = 0, the vertical deflection is 0. As t approaches infinity, the vertical deflection is equal to the vertical deflection calculated by the traditional elastic half-space theory. The viscosity coefficient η generally ranges from 0 to 1, which limits the overall vertical deformation. It also causes the theoretically calculated peak value of deflection time-history to reduce and the phase to extend to the right.


### 3.3. Establishment of Modulus Back-Calculation Procedure

As shown in equation (5), the back-calculation parameters are modulus value E and the viscosity coefficient η. Due to the complexity of equation (5), it is difficult to obtain the analytic solution to determine the back-calculation parameters. Therefore, the least square method was used in this study, namely the optimization algorithm iterated by(6)minΔ=K∑i=1nWti−wti,β2,where K is the accuracy control parameter, W is the measured value of vertical deflection (mm), β=E,η is the back-calculation parameter vector, and i and n are the data number and the total data number of deflection time-history measured by PFWD, respectively.The modified Gauss–Newton method was used in the iterative process, that is, whenminΔ=fβ, equation (6) was transformed into equation (7) as follows:(7)fβ=KrTβrβ,where rβ is the residual column vector of measured and theoretical values.In the calculation, the initial calculation parameter vectorβ0 was first set. An initial modulus and a viscous coefficient value are generally set by experience. The next iteration parameter vector was determined by(8)βk+1=βk+λkAk,where k is the iteration number, k = 0, 1, 2, … , and λk is the step size, which can be 0.01∼0.05 by experience.In order to improve the iteration speed, the optimal value of step size can also be searched byminβk+λkAk, where Ak calculated by Ak=−JrTJr−1JrTrβk is the direction vector and Jr is the Jacobi matrix when β=βk.The iterative cycle process was controlled by(9)βk+1−βkT⋅10≤ε,where T is the vector transpose and ε is the iterative control parameter.The final iteration valueβg was calculated by equation (10), that is, the optimal back-calculation value of modulus and viscosity coefficient.(10)βg=βk+1+βk2.Because the least-square optimization algorithm is based on the back-calculation result of global extremum in the function definition domain, the results may be quite different from the actual situation. Therefore, when the back-calculation program was compiled, the specific conditions of high plasticity soil modulus detected by PFWD were dealt with as follows. Since the data collection period is 0∼60 ms, the PFWD will vibrate on the rubber block, thus driving the bearing plate to have a small jump. Therefore, the 0∼25 ms peak waveform data were taken. The calculation point wasr = 0 m, that is, the center of load acting surface. The accuracy control parameter K is 0.5, and the iterative control coefficient ε is 0.01. These two control values can be adjusted according to the accuracy requirements.In addition, equation (8) was revised into equation (11) to control the weight of the back-calculation modulus. When the phase of calculated deflection to time-history is greater than the measured deflection time-history, or the calculated deflection peak value is less than the measured deflection value, the iteration of parameter η will stop and the iteration of parameter E will continue because the iteration fails to meet the requirements of equation (9).(11)βk+1=βk+λkAkB,where B is the iterative weight control vector calculated by (12)B=1,1,Eβk≤γ,1,0,Eβk>γ,where Eβk is the limit condition of peak error of measured deflection and calculated deflection (Eβk=Wmax−wmaxβk), Wmax is the measured deflection peak (mm), wmaxβk is the calculated deflection peak (mm), and γ is the deflection peak error control parameter (mm), which can be adjusted to change the accuracy of the final back-calculation according to the users’ requirement. In this study, γ=0.002mm.After obtaining the back-calculation solution, the deflection was calculated and the calculation accuracy was evaluated by(13)R2=1−∑i=1nWti−wti,βk2∑i=1nWti2−∑i=1nWti2/n,where R2 is the accuracy factor.According to the method proposed in this study, the back-calculation program of soil modulus considering the effect of deformation hysteretic effect under impact load was compiled by MATLAB. It applies to both FWD and PFWD back-calculations. The proposed back-calculation process is shown in Figure8.Figure 8
Proposed back-calculation method.


## 3.1. Selected Viscoelastic Theoretical Model

It is well known that under a sinusoidal dynamic load, the strain of viscoelastic material changes periodically, and the strain response lags behind the stress response. The viscoelastic theoretical model has been commonly applied to dynamic response analysis of asphalt pavement under the traffic load. The deformation hysteresis effect of a high plasticity soil embankment under impact load is similar to the mechanical behavior of viscoelastic materials mentioned earlier. Therefore, based on the viscoelastic constitutive model, a high plasticity soil modulus back-calculation method was established in this research.The Kelvin and Maxwell viscoelastic models are relatively simple, whereas the three-parameter solid model, Burgers model, and generalized Maxwell model are more complex. The Maxwell model is a series of spring-dampers that are primarily used to simulate the stress relaxation process of viscoelastic materials. The three-parameter solid model, Burgers model, and a generalized Maxwell model are constructed by a spring-damper parallel and series. The calculated parameters were incorporated to better describe the viscoelastic and residual deformation properties of material under dynamic load. However, the multiparameter model greatly increases the difficulty and iteration time of the back-calculation, but the results may not converge. In the PFWD test, the loading and unloading time of drop hammer impact is very short so that the plastic deformation cannot occur in time. The vertical deformation measured by PFWD is mainly rebound deformation. However, only the recoverable deformation was used to determine the elastic modulus, and the permanent deformation part was eliminated by the repeated dropping of the hammer. Therefore, it is not necessary to consider the residual deformation or adopt the complicated multiparameter viscoelastic model to calculate the high plasticity soil modulus.In this study, the Kelvin model was selected, and its concise structure was conducive to an efficient iteration of subsequent MATLAB program. The mechanical response of subgrade can be described by this model to meet the accuracy requirements.The Kelvin model, a viscoelastic model, was composed of a spring-damper in parallel. It is simple in form and easy to iterate efficiently. Its mathematical description is as follows:(2)σ=Eε+ηdεdt,where σ is the stress, ε is the strain, E is the elastic modulus, and η is the viscosity coefficient.The load-time curve obtained from the PFWD impact load can be considered as a semisine curve. Assuming that the mathematical expression of curve isσ¯sinωt, equation (3) was obtained by substituting σ¯sinωt into equation (2) and ignoring the rapid attenuation index items associated with time t:(3)σσ¯2−2cosδσσ¯εε¯+εε¯2=sin2δ,where σ¯ is the stress amplitude, ε¯ is the strain amplitude, δ is the phase difference, and tanδ=ηω/E.According to equation (3), the stress-strain curve under the action of loading and unloading conditions can be obtained, as shown in Figure 7. It could be seen that when the viscosity coefficient was introduced, the loading stress first reached the maximum value, and the corresponding strain did not reach the maximum. However, the maximum strain was observed during the unloading process due to the hysteresis. The stress-strain curve is very similar to the PFWD-measured load-deflection curve, as shown in Figure 4(b). Therefore, the Kelvin model was adopted to analyze the dynamic response of high plasticity soil under impact load.Figure 7
Stress-strain curve of Kelvin model.


## 3.2. Quasi-Static Back-Calculation Method

The dynamic and quasi-static methods were used for simulating dynamic response analysis of elastic systems under dynamic loading.In the dynamic back-calculation [20], combined the motion differential equation of axisymmetric space problem in the cylindrical coordinate system with the physical equation of stress expressed by displacement, the Navier dynamic equation can be obtained. The viscoelastic theory model was introduced into the Navier dynamic equation, and the explicit solution of vertical deflection was derived by Laplace and Hankel back-calculation. The complex mixed-boundary-value problem was also considered. Since the structure of implicit partial differential equations was complicated in back-calculation, time-consuming numerical algorithms were adopted to simplify it.For this research, the quasi-static back-calculation method is to take the explicit solution of vertical deflection on the surface of the elastic half-space body under the static load as Laplace transformation and to replace the elastic parameters in the equation with viscoelastic operators. In this method, the dynamic load data were used to calculate the static parameters. Subsequent back-calculation results verified that the quasi-static method met the requirements for back-calculation of high plasticity soil embankments. Therefore, the modulus calculated by this method and the PFWD dynamic load data described in this study is static modulus, which is consistent with the static modulus measured by the PLT.The vertical deflection on the surface of elastic half-space body under the action of an axial symmetric load was calculated as follows:(4)wr=21−μ2E∫0∞J0ξrp¯ξdξ,where J0(x) is the zero-order Bessel function, p¯ξ is the zero-order Hankel integral transformation of the load for PFWD rigid bearing plate, p¯ξ=pRsinξR/2ξ; p is the load plate pressure (kN⋅m−2), R is the radius of the load plate (m), and r is the radial distance from the load center of calculation point (m).After the Laplace transformation of equation (4) was performed, the viscoelastic operator was used to replace the elastic modulus, and the PFWD load function p¯ξ was substituted into the equation for inverse Laplace transformation. Since the Poisson ratio, which changes less, has little influence on the calculation result, it was not considered in the calculations. The vertical explicit deflection equation of viscoelastic half-space surface under the action of the PFWD rigid loading plate was obtained as follows:(5)wr,t=L−1w¯r,s=1−μ2Fπ⋅Es⋅R1−e−E/ηt∫0∞J0ξrsinξRξdξ,where E(s) = E + ηs is the viscoelastic operator; F is the concentrated force generated by the PFWD (KN) acting on the surface of the measuring point, which is obtained by the load time-history curve; R and r have the same meaning as equation (4); and s is the Laplace transformation parameter of time t.In equation (5), it is shown that when t = 0, the vertical deflection is 0. As t approaches infinity, the vertical deflection is equal to the vertical deflection calculated by the traditional elastic half-space theory. The viscosity coefficient η generally ranges from 0 to 1, which limits the overall vertical deformation. It also causes the theoretically calculated peak value of deflection time-history to reduce and the phase to extend to the right.


## 3.3. Establishment of Modulus Back-Calculation Procedure

As shown in equation (5), the back-calculation parameters are modulus value E and the viscosity coefficient η. Due to the complexity of equation (5), it is difficult to obtain the analytic solution to determine the back-calculation parameters. Therefore, the least square method was used in this study, namely the optimization algorithm iterated by(6)minΔ=K∑i=1nWti−wti,β2,where K is the accuracy control parameter, W is the measured value of vertical deflection (mm), β=E,η is the back-calculation parameter vector, and i and n are the data number and the total data number of deflection time-history measured by PFWD, respectively.The modified Gauss–Newton method was used in the iterative process, that is, whenminΔ=fβ, equation (6) was transformed into equation (7) as follows:(7)fβ=KrTβrβ,where rβ is the residual column vector of measured and theoretical values.In the calculation, the initial calculation parameter vectorβ0 was first set. An initial modulus and a viscous coefficient value are generally set by experience. The next iteration parameter vector was determined by(8)βk+1=βk+λkAk,where k is the iteration number, k = 0, 1, 2, … , and λk is the step size, which can be 0.01∼0.05 by experience.In order to improve the iteration speed, the optimal value of step size can also be searched byminβk+λkAk, where Ak calculated by Ak=−JrTJr−1JrTrβk is the direction vector and Jr is the Jacobi matrix when β=βk.The iterative cycle process was controlled by(9)βk+1−βkT⋅10≤ε,where T is the vector transpose and ε is the iterative control parameter.The final iteration valueβg was calculated by equation (10), that is, the optimal back-calculation value of modulus and viscosity coefficient.(10)βg=βk+1+βk2.Because the least-square optimization algorithm is based on the back-calculation result of global extremum in the function definition domain, the results may be quite different from the actual situation. Therefore, when the back-calculation program was compiled, the specific conditions of high plasticity soil modulus detected by PFWD were dealt with as follows. Since the data collection period is 0∼60 ms, the PFWD will vibrate on the rubber block, thus driving the bearing plate to have a small jump. Therefore, the 0∼25 ms peak waveform data were taken. The calculation point wasr = 0 m, that is, the center of load acting surface. The accuracy control parameter K is 0.5, and the iterative control coefficient ε is 0.01. These two control values can be adjusted according to the accuracy requirements.In addition, equation (8) was revised into equation (11) to control the weight of the back-calculation modulus. When the phase of calculated deflection to time-history is greater than the measured deflection time-history, or the calculated deflection peak value is less than the measured deflection value, the iteration of parameter η will stop and the iteration of parameter E will continue because the iteration fails to meet the requirements of equation (9).(11)βk+1=βk+λkAkB,where B is the iterative weight control vector calculated by (12)B=1,1,Eβk≤γ,1,0,Eβk>γ,where Eβk is the limit condition of peak error of measured deflection and calculated deflection (Eβk=Wmax−wmaxβk), Wmax is the measured deflection peak (mm), wmaxβk is the calculated deflection peak (mm), and γ is the deflection peak error control parameter (mm), which can be adjusted to change the accuracy of the final back-calculation according to the users’ requirement. In this study, γ=0.002mm.After obtaining the back-calculation solution, the deflection was calculated and the calculation accuracy was evaluated by(13)R2=1−∑i=1nWti−wti,βk2∑i=1nWti2−∑i=1nWti2/n,where R2 is the accuracy factor.According to the method proposed in this study, the back-calculation program of soil modulus considering the effect of deformation hysteretic effect under impact load was compiled by MATLAB. It applies to both FWD and PFWD back-calculations. The proposed back-calculation process is shown in Figure8.Figure 8
Proposed back-calculation method.


## 4. Verification of the Proposed Modulus Back-Calculation Method

The load-deflection time-history data obtained from L1-10 and U1-10 were used to back-calculate the elastic modulus with the VEM. The LEM and VEM were used to obtain the back-calculated modulus for these locations which was further compared with that measured by the PLT. The moduli of L4 and U4 obtained from different methods are also listed in Table1. The estimated moduli for the untreated (lower embankment) and cement-treated (upper embankment) soils are 17.2 MPa and 192.3 MPa, respectively. The proposed VEM method overestimated the moduli of untreated and treated soils by 13.9% and 3.8%, respectively. However, the conventional LEM overestimated the moduli of untreated and treated soils by 67.5% and 15.3%, respectively. The VEM can greatly reduce the error in estimating the modulus as it can accurately predict the time required to reach the peak deflection obtained from the PFWD. Compared with the actual in situ data, the moduli obtained from the VEM and LEM were used to regenerate the deflection-time-history data.Table 1
Comparison of the LEM, VEM, and measured results.
Measuring point
L4
U4
LEM
VEM
Measured
LEM
VEM
Measured
Modulus (MPa)
25.3
17.2
15.1a
213.7
192.3
185.3a
R
2
0.693
0.935
NA
0.826
0.945
NA
Load
Peak (kN)
NA
NA
7.313b
NA
NA
7.160b
Time (ms)
NA
NA
6b
NA
NA
8b
Deflection
Peak (μm)
845
872
877b
98
99
99b
Time (ms)
6
9
9b
8
9
9b
Note:R2 is the accuracy coefficient, which is calculated by equation (13); a is the test result obtained from PLT; b is the test result obtained from PFWD.In Figure9, a comparison of the time-history curve of L4 is shown. It can be seen that the deflection time-history calculated by the VEM is in good agreement with the measured deflection time-history. According to equation (13), the calculation accuracy coefficients measured by the VEM and LEM are 0.935 and 0.693, respectively.Figure 9
Comparison of measured and back-calculated deflection-time curves of L4.In Figure10, a comparison of the time-history curve of U4 is shown. It can be seen that the deflection time-history calculated by the VEM and LEM are in good agreement with the measured deflection time-history. For the cement-treated soil, the calculation accuracy coefficients measured by the VEM and LEM are 0.945 and 0.826, respectively. Under the same loading condition, for the stiffer soil, the vertical deflection is smaller, which reduces the effect of filling viscosity. Therefore, the moduli obtained from LEM and VEM are very close to each other.Figure 10
Comparison of measured and back-calculated deflection-time curves of U4.Back-calculated moduli of all 10 test points on the top of lower embankment performed by the VEM are shown in Figure11(a). The maximum, minimum, and average back-calculation modulus values of 10 test points are 21.5 MPa, 17.2 MPa, and 19.5 MPa, respectively. The modulus obtained from the VEM is closer than that obtained from the LEM. The average absolute error is only 2.5 MPa, and the average relative error is14.8%. It has been proven that the quasi-dynamic analysis method can be used to calculate the static parameters according to the dynamic response of embankments under impact load (load-deflection time-history data). The static modulus measured by VEM is closer to that measured by the PLT.Figure 11
Comparison of back-calculated moduli using VEM and measured moduli using PLT. (a) Results of L1–L10. (b) Results of U1–U10.
(a)
(b)In Figure11(b), the measured and back-calculated results of all 10 test points on the top of the K110 + 180∼K110 + 280 trial section of upper embankment are shown. The statistical results show that the maximum, minimum, and average modulus values of 10 test points back-calculated by the VEM are 220.6 MPa, 191.6 MPa, and 205.8 MPa, respectively. The average absolute error is 8.5 MPa, and the average error is only 4.3%. It is shown that when the soil modulus is high, the average relative errors of moduli back-calculated by the VEM and LEM are both reduced, but the former is smaller.According to the field test results of trial section of high plasticity soil embankment, the relative errors of back-calculated modulus by the VEM and measured modulus are both within 15% in the lower and upper embankments with low stiffness and stiffer soils, respectively. Therefore, the back-calculation modulus measured by the VEM can be used as the parameter for pavement design.


## 5. Case Study

The soil modulus of upper embankment obtained from the VEM is statistically lower than that obtained from the LEM. Ignoring the viscoelastic nature of embankment may lead to the overestimation of stiffness, which has a significant influence on the design life of a pavement. A hypothetical case study was performed to analyze the effect of soil modulus on the design life of a flexible pavement consisting of an asphalt layer, base layer, and subgrade (Figure12). Rutting, fatigue, and low-temperature cracking are the principal types of diseases that should be considered for the design of flexible pavement [30]. Typical modulus values obtained from the literature were used for the asphalt and base layers. Soil moduli of upper embankment obtained from both the LEM and VEM were used as subgrade moduli for comparison.Figure 12
Pavement section. Note: EA is the elastic modulus of the asphalt layer; EB is the elastic modulus of the base layer; EL is the average modulus (U1–U10) back-calculated by LEM; EV is the average modulus (U1–U10) back-calculated by VEM.Considering the soil moduli obtained from the LEM and VEM, pavement analyses were performed with FPS 21 software to compare the cracking and rutting life of flexible pavement. A low-volume road expected to carry 0.8 million equivalent standard axle loads (ESALs) throughout its design life of 20 years was considered. Design data used for the pavement analysis are listed in Table2.Table 2
Design data.
Parameters
Value
Parameters
Value
Length of analysis
20
Initial ADT (vehicles/day)
800
Design confidence level (95%)
C
Final ADT (vehicles/day)
1200
Initial serviceability index
4.5
Total 80 kN ESALs (millions)
0.8
Final serviceability index
2.5
Percentage of truck in ADT
7%


### 5.1. Fatigue Cracking Criteria

The fatigue cracking of flexible pavement can be determined by the following equation:(14)Nf=f1εt−f2EA−f3,where Nf is the allowable number of load repetitions; εt is the horizontal tensile strain acting at the bottom of asphalt layer; EA is the elastic modulus of the asphalt layer; and f1, f2, and f3 are constant values proposed by the agencies that are listed in Table 3.Table 3
Coefficients for determining the cracking life.
Agency
f
1
f
2
f
3
Asphalt Institute
0.0796
3.291
0.854
Shell equation
0.0685
5.671
2.363
Transport and road research laboratory (TRRL)
1.66 × 10−10
4.320
0Cracking life determined in terms of traffic volume depends on the thickness of base layer, as shown in Figure13.Figure 13
Fatigue cracking life of flexible pavement. (a) Asphalt Institute method. (b) Shell method. (c) TRRL method.
(a)
(b)
(c)A significant reduction in the design life was observed if the soil modulus of upper embankment obtained from the VEM were used instead of that obtained from the LEM. In this case study, the soil moduli obtained from the LEM and VEM are 222 MPa and 206 MPa, respectively. The Asphalt Institute, Shell, and TRRL methods were used for the analysis, which showed that with a reduction of 8% in the soil modulus, the cracking life measure by the aforementioned three methods decreased by 12.6%, 20.8%, and 15.5%, respectively.


### 5.2. Rutting Criteria

Rutting, a common phenomenon of flexible pavement, is indicated by the permanent deformation or rut depth that occurs along the wheel paths. Design life (Nd: allowable load repetition) can be calculated by the following equation:(15)Nd=f4εv−f5,where εv is the vertical compressive strain acting on the top of subgrade layer and f4 and f5 are constant values proposed by different agencies that are listed in Table 4.Table 4
Coefficients for determining the rutting life.
Agency
f
4
f
5
Asphalt Institute
1.36 × 109
4.477
Shell (95% reliability)
1.05 × 10−7
4.000
UK transport and road research (UTRR)
6.18 × 10−8
3.950The rutting life obtained from these methods is shown in Figure14.Figure 14
Rutting life of flexible pavement. (a) Asphalt Institute method. (b) Shell method. (c) UTRR method.
(a)
(b)
(c)The compressive strain acting on the top of the subgrade layer is directly related to the soil modulus of upper embankment and the thickness of the base layer. Reducing the soil modulus led to an increase in the vertical strain, which eventually reduced the pavement rutting life. With a reduction of 8% in the soil modulus, the pavement rutting life decreased by 25% (base thickness = 15 cm). Based on the Asphalt Institute’s method, when the Shell and UTRR methods were used, the rutting life decreased by 18.4% and 19.6%, respectively.


## 5.1. Fatigue Cracking Criteria

The fatigue cracking of flexible pavement can be determined by the following equation:(14)Nf=f1εt−f2EA−f3,where Nf is the allowable number of load repetitions; εt is the horizontal tensile strain acting at the bottom of asphalt layer; EA is the elastic modulus of the asphalt layer; and f1, f2, and f3 are constant values proposed by the agencies that are listed in Table 3.Table 3
Coefficients for determining the cracking life.
Agency
f
1
f
2
f
3
Asphalt Institute
0.0796
3.291
0.854
Shell equation
0.0685
5.671
2.363
Transport and road research laboratory (TRRL)
1.66 × 10−10
4.320
0Cracking life determined in terms of traffic volume depends on the thickness of base layer, as shown in Figure13.Figure 13
Fatigue cracking life of flexible pavement. (a) Asphalt Institute method. (b) Shell method. (c) TRRL method.
(a)
(b)
(c)A significant reduction in the design life was observed if the soil modulus of upper embankment obtained from the VEM were used instead of that obtained from the LEM. In this case study, the soil moduli obtained from the LEM and VEM are 222 MPa and 206 MPa, respectively. The Asphalt Institute, Shell, and TRRL methods were used for the analysis, which showed that with a reduction of 8% in the soil modulus, the cracking life measure by the aforementioned three methods decreased by 12.6%, 20.8%, and 15.5%, respectively.


## 5.2. Rutting Criteria

Rutting, a common phenomenon of flexible pavement, is indicated by the permanent deformation or rut depth that occurs along the wheel paths. Design life (Nd: allowable load repetition) can be calculated by the following equation:(15)Nd=f4εv−f5,where εv is the vertical compressive strain acting on the top of subgrade layer and f4 and f5 are constant values proposed by different agencies that are listed in Table 4.Table 4
Coefficients for determining the rutting life.
Agency
f
4
f
5
Asphalt Institute
1.36 × 109
4.477
Shell (95% reliability)
1.05 × 10−7
4.000
UK transport and road research (UTRR)
6.18 × 10−8
3.950The rutting life obtained from these methods is shown in Figure14.Figure 14
Rutting life of flexible pavement. (a) Asphalt Institute method. (b) Shell method. (c) UTRR method.
(a)
(b)
(c)The compressive strain acting on the top of the subgrade layer is directly related to the soil modulus of upper embankment and the thickness of the base layer. Reducing the soil modulus led to an increase in the vertical strain, which eventually reduced the pavement rutting life. With a reduction of 8% in the soil modulus, the pavement rutting life decreased by 25% (base thickness = 15 cm). Based on the Asphalt Institute’s method, when the Shell and UTRR methods were used, the rutting life decreased by 18.4% and 19.6%, respectively.


## 6. Conclusions

In this study, based on the viscoelastic model, a modulus back-calculation method was proposed to estimate the static soil modulus which used the load-deflection/time-history data under the impact load. The effectiveness of the VEM was verified by the in situ PFWD and PLT performed on low and high stiffness soil layers. The effect of moduli obtained from the LEM and VEM on the pavement performance was investigated through a hypothetical case study. The following conclusions can be drawn:(1)
Under an impact load, the peak value of deflections lags behind the peak value of load, and the viscoelastic hysteresis of load-deflection curve is noticeable. The subgrade has significant deformation hysteresis and load-displacement nonlinearity.(2)
The nonlinear dynamic response characteristic of soil was ignored by the LEM, and the modulus was overestimated. For a stiffer soil (>200 MPa), the average error is around 12%, but for low stiffness soil (<25 MPa), the average error is more than 50%. The average error of LEM decreases with the increase in soil modulus.(3)
The loading frequency for the PLT is 0 Hz. However, a frequency of about 4 Hz was used by the LEM, which led to the overestimation of modulus. This frequency has certain influence on the measurement and back-calculation of soil modulus.(4)
The VEM of PFWD modulus back-calculation can reflect the nonlinear dynamic response characteristic of soils. It can also calculate the static modulus by incorporating the dynamic load data. For a stiffer soil, the average error is less than 5%; for the low stiffness soil, the average error is less than 15%.(5)
Pavement analysis in the hypothetical case study shows that compared with the modulus obtained from the VEM, the modulus obtained from the LEM overestimates the fatigue-cracking and rutting life of flexible pavement.



---

*Source: 1316341-2019-11-06.xml*