INTRODUCTION
The modulus of elasticity is the major parameter in the remaining life analysis
and in the design of new pavements and overlays. This critical value is computed
through the backcalculation method, which uses the FWD raw data (Tawfiq,
2003). Most of the common backcalculation algorithms, which use the peak
values of load and deflection at each sensor during every drop, assume that
the FWD load is applied statically. The pavement system is modeled as a layered
elastic system with linear or nonlinear (stressdependent) materials. A forward
analysis subroutine is used to calculate the theoretical deflections under the
known load. Attempts are made to converge these calculated deflections to the
measured ones. In other words, the moduli are predicted by minimizing the error
between the measured and calculated deflections in every step of iteration.
But in reality, the FWD applies an impulse load to the pavement that generates
body waves and surface waves, which travel at finite velocities and are recorded
at different times by the geophones (Lytton, 1989). This
full time history data of both load and deflections can also be measured during
the FWD test. This additional data can be used to backcalculating other parameters
of pavement materials.
BACKCALCULATION OF LAYER MODULI, STATIC AND DYNAMIC FEATURES
Classical pavementsubgrade system analysis is based on multilayer linear
elastic analysis, infinite pavement layer dimensions and a semiinfinite subgrade.
Linear elastic response assumptions are valid when there is no discontinuity,
for example at joints in rigid pavements, or cracks. The assumed static loading
conditions are also a simplification of reality (William, 1999).
The duration of loading runs from 2065 msec depending on the type of FWD device.
The FWD provides excellent structural information if it is fully utilized. From
a review of the sensor peak values for each measurement location, it becomes
clear that points farther from the center of load plate attain their peak values
later than locations closer to the plate. This time differences known as phase
difference of deflections and is a measure of the velocity of propagation of
the shockwaves (Fig. 1). These types of time history data
contain useful information that if well utilized would contribute to improved
accuracy of the structural pavement evaluation. The use of only peak values
discards this potentially very useful information (Matsui
and Hachia, 2006).
Static Feature
Objective Function
The objective of most backcalculation programs is to determine a set of
moduli that will minimize an error term between the computed deflection and
the measure deflection (Chua, 1989). Because random errors
are unavoidable, the number of measuring deflection sensors must be greater
than the number of unknowns (for example, moduli in the linear elastic case
or material parameters for the nonlinear case) to be backcalculated (Uzan
et al., 1989).
Since, the accuracy of the sensors, which is a major source of error, is expressed
in relative terms, the objective function to minimize should therefore be expressed
in relative terms. The objective function is (Uzan et al.,
1989):
ε^{2} 
= 
Squared error 
w_{i}^{m} 
= 
Measured deflection at sensor i 
w_{i}^{c} 
= 
Computed deflection at sensor i 
S 
= 
The number of sensors 
Backcalculation Methods
The computed deflection w_{i}^{c} can be expressed as follows
(Uzan et al., 1989):
w_{i}^{c} = f_{i}(X_{j}) 
(2) 
X_{j} is the unknown variables, including modulus of elasticity,
poisson’s ratio, thickness of layer, pressure, contact area,…
(j varies from 1 to n).
Generally, there are two groups of backcalculation programs. The first is the
iterative approach in which a forward calculation scheme (typically a linear
elastic program) is used within the iterative process (Chua,
1989). In other words, any solution to Eq. 1 calls for
a solution to Eq. 2. The number of calls depends on the minimization
algorithm used. So, the pattern search technique requires a hundreds of calls
of the deflection computation program. A simplified description of the iterative
process used for adjusting the modulus values is shown in Fig.
2. This illustration is for one deflection and one layer. For multiple layers
and deflections, the solution is obtained by developing a set of equations that
define the slope and intercept for each deflection and each unknown layer modulus
as follows (Van Cauwelaert et al., 1989; Tawfiq,
2003):
Log (deflection_{j}) = A_{ij}+
S_{ji} (log E_{i}) 
(3) 
A 
= 
Intercept 
S 
= 
Slope 
j 
= 
1, 2,…, ND (ND = No. of deflections) 
i 
= 
1, 2,…, NL (NL = No. of layers with unknown moduli) 
The data base method uses forward calculation scheme to build a data base from
which regression equations are either formulated to determine the layer moduli
or are used within interpolation techniques to compute the deflections, thus
avoiding the use of a forward calculation scheme in the iterative process (Uzan
et al., 1989).
In the case of linear elasticity and a circular contact area, Eq.
2 can be written as (Uzan et al., 1989):
p 
= 
Pressure 
a 
= 
Radius of contact area 
E_{i} 
= 
Modulus of elasticity of the ith layer on the subgrade 
E_{sg} 
= 
Subgrade modulus of elasticity 
This equation is similar to the well known equation for the deflection of two
layer systems, as shown in Fig. 3 (Tawfiq,
2003):
• 
for flexible pavement 
• 
for rigid pavement 
F_{2} is a function of E_{1}/E_{2} and h_{1}/a
. The value of f_{i}, for sensor i, which is a function of modular ratios,
can be obtained from the layered system program and used as a data base (Tawfiq,
2003).
Backcalculation Programs
Most of the common backcalculation programs with elastostatic basis
use an elastic layer program for forward analysis and an error optimization
algorithm. In this study, the programs MODULUS 6.0, ELMOD 5.0 and EVERCALC
5.0 have been used to perform traditional back analysis.
Description of MODULUS6.0
This program was developed at Texas Transportation Institute (TTI) for the
Texas DOT and uses the database method for backanalysis. It uses WESLEA program
as a forward calculation subroutine. WESLEA is based on the multilayer linear
elastostatic theory that is traditionally used for the purposes of flexible
pavement analysis (William, 1999). The program uses WESLEA
to generate a data base of deflection bowls by assuming different modular ratios.
A pattern search technique is then used to determine the set of layers moduli
that produce a deflection basin that fits the measured one (William,
1999).
The problem with the program is that the units are in inlb units only
and that it only handles 7 sensors.
Up to four unknown layers, with or without bedrock are allowed. It is enough
for the user to define the material type and the thickness for every layer.
The program suggests the moduli ranges and poisson’s ratios and also the
depth to bedrock and includes them in the analysis. The subgrade depth can be
altered by the user if site specific information is available (William,
1999).
Description of ELMOD5.0
ELMOD is an acronym for Evaluation of Layer Moduli and Overlay Design. This
program was developed by Dynatest Consulting Inc. and accepts DYNATESTFWD format.
The theory of OdemarkBoussinesqe transformed section is the basis of forward
analysis. Besides, it uses an iterative method for backcalculation (William,
1999). There are 3 modes of backcalculations available (Elmod
5 Quick Start Manual, 2006):
• 
Radius of curvature: The outer geophone readings are used
to determine nonlinear characteristics of subgrade and the inner
geophones to determine the upper pavement layer moduli. The stiffness
of remaining layers is then calculated based on the overall pavement
response to the applied load. This approach analysis up to 4 layers 
• 
Deflection basin fit: The methodology starts with a set of
estimated moduli for the pavement structure. The theoretical deflection
bowl for this pavement structure is calculated. The error between
the measured and calculated deflections is then assessed. The moduli
in the structure are increased/decreased by a small amount (typical
10%) and if the error in either of these bowls is less than the original
bowl, this is taken to be a better solution. Up to 5 layers can be
estimated by this mode 
• 
FEM/LET/MET: With this option backcalculation may be carried out
either with Finite Element Method (FEM), Linear Elastic Theory (LET) or
the Method of Equivalent thickness. FEM makes use of a modified version
of an axial symmetric finite element program, originally developed by Wilson
at University of California by Duncan et al. (1968).
LET makes use of WESLEA and MET is similar to the method used in the option
Basin Fit, but with a simpler use of adjustment factors 
Temperature correction is also engaged. Moduli range and poisson’s ratio
(equal to 0.35 for all layers in every case) is determined by the software internally
(Elmod 5 Quick Start Manual, 2006).
The major assumption is that the layers are homogeneous, isotropic and linear
elastic, except the subgrade which can have a nonlinear definition (Elmod
5 Quick Start Manual, 2006):
^{} 
(5) 
σ_{1} 
= 
Major principal stress from external loading 
p_{a} 
= 
Reference stress, often taken equal to atmospheric pressure (0.1
MPa) 
C and n are constants. n varies from 0 (linear elastic material) to 0.5.
Description of EVERCALC5.0
EVERCALC5.0 was developed by Dr. Joe Mahoney at the University of Washington
for Washington DOT. This program uses the WESLEA computer code for the forward
analysis and a modified GaussNewton algorithm for solution optimization. It
can handle up to 5 layers, 10 sensors and 12 drops per station (William,
1999).
The program terminates when one or more of the following conditions are satisfied
(Everseries User’s Guide, 2005):
w_{i}^{m} is the measured deflection and w_{i}^{c}
is the calculated deflection at sensor i. S is the number of sensors and
n is the number of layers
E_{ki} and E_{(k+1)i} are the ith layer moduli at the
kth and (k+1)th iteration and m is the number of layers with unknown moduli

Number of iterations has reached the maximum number of iterations.
At every iteration a maximum of (m+1) calls to WESLEA is made, where
m is the number of layers with unknown moduli 
Dynamic Feature
System Identification Methods
System identification methods were developed by electrical engineers who
were interested in determining the characteristics of a filter by using the
input and output signals of the filter and an assumed model of the filter. The
characteristics of the filter are changed systematically using a search technique
until the model produces an output that is acceptably close to that of the filter.
The procedure is in fact exactly analogous to what is being done in backcalculating
the moduli of pavements (Lytton, 1989).
An impulse load is applied to the pavement at a remote distance and the surface
motion is sensed as signals received by two sensors, which are separated by
a distance, x. The pavement between the two sensors is regarded as the filter,
the characteristics of which are to be determined. The unknowns in each layer
are its modulus, percent damping, thickness and unit weight (Lytton,
1989).
Description of DBSID (Dynamic Backcalculation with System Identification
Method): The Fourier transform is commonly used and widely implemented as
an integral transform in signal processing which decomposed e.g., a function
of time into a series of harmonics. There are several efficient algorithms for
implementing the discrete Fourier transform, commonly referred to as Fast Fourier
Transform (FFT) (Westover and Guzina, 2006). The FWD load
is first decomposed into its frequency components using the FFT. Then the transfer
function, which defines the response of the pavement system to a steady state
unit load, is evaluated. Finally, the transfer function is multiplied by the
FFT of the load to determine the Fourier transform of the displacements (Fernando
and Liu, 2002).
The dynamic analysis presented by DBSID is based on the forward model implemented
in the FWDDYN program. Based on research by Uzan et al.
(1989), the developers of DBSID followed a timedomain fitting approach
in the dynamic backcalculation. This approach was readily implemented using
the FWDDYN program, which predicts the displacement history for each FWD sensor,
given the material properties and thickness of each pavement layer. FWDDYN
uses Fourier superposition to predict pavement response due to the impulse load
by the FWD. In this method, FWDDYN performs an inverse FFT on this Fourier
transform to determine the time history of the displacements for each FWD sensor.
The developers of DBSID, coupled a system identification routine to the FWDDYN
program. The Asphalt Concrete (AC) layer is modeled as a viscoelastic (or damped
elastic) material and the base and subgrade layers are modeled as damped elastic
materials (Fernando and Liu, 2002).
Here, point minimizing of the errors between discrete values of deflections
is not the concern, but the best fit of the deflection histories.
Pavement Dynamic Material Model
Viscoelastic Material Model (User’s Manual for DBSID, 1993)
Asphalt Concrete (AC) material is a matrix of solid aggregate particles
with asphalt binder. The bonding between the asphalt binder and the aggregate
particles results in viscoelastic behavior of AC mixture. In the DBSID
procedure, the viscoelastic behavior is modeled using the following threeparameter
model that is considered realistic for loading times representative of
highway traffic speed:
D(t) = D_{0}+ D_{1} t^{m} 
(8) 
D(t) 
= 
The creep compliance function, defined as D(t) = e (t)/s_{0} 
e(t) 
= 
The total axial strain measured at time t 
s_{0} 
= 
The applied constant stress in a creep compliance test 
D_{0}, D_{1} and m 
= 
Model parameters 
When t = 0, D(0) = D_{0}, from the threeparameter viscoelastic
model. This means that D_{0} (D_{0} = 1/E_{0})
stands for the elastic response at a very small loading time. For the
asphalt concrete layer, the following parameters are used in backcalculation:
H 
= 
Thickness of asphalt concrete layer 
g 
= 
Unit weight of asphalt concrete 
v 
= 
Poisson’s ratio 
D_{0} 
= 
1/E_{0}, representing the elastic response due to the solid
matrix 
D_{1} 
= 
Creep compliance constant representing the nonlinear viscous response 
m 
= 
Exponent for nonlinear time dependence 
DampedElastic Material (User’s Manual for DBSID, 1993)
For pavement layers that are modeled as dampelastic, only two parameters
are needed to define the response of the material: the elastic modulus
(E) and the damping ratio (ξ). The damping ratio is input to the
program and is not backcalculated. Only the elastic modulus is backcalculated.
If the surface layer is modeled as viscoelastic, the values D_{0},
D_{1} and m are backcalculated. In the other damped elastic layers
the moduli are backcalculated. So, together with the depth of subgrade,
6 unknowns will be defined through backcalculation.
FIELD DATA COLLECTION
The FWD data used in the analysis was collected on the freeways and airport
site listed below:
• 
Eivanekey Garmsar, km 00.00 to 26.55 
• 
Garmsar Semnan, km 00.00 to 104.13 
• 
ZanjanTabriz, km 205.00 to 236.00 
• 
Rafsanjan Airport, in three different parts: Taxiway, Runway and
Apron 
It is noteworthy that in the last item, the HWD test has been performed.
COMPARISON BASIS
In order to evaluate the results, a comparison basis is essential. The subgrade
layers are the most complex layers in the pavement structure due to their physical
nature and construction practices. The backcalculation process is based on the
assumption that surface deflection at a certain offset is characteristic of
the elastic modulus at a certain depth (William, 1999).
Several of the backcalculation programs start their analysis from the outer
geophone which determines the resilient modulus of subgrade. So, accurate assumption
of subgrade is vital to the rest of the analysis.
In this study, the results of soil mechanic tests for subgrade layers
have been used as a groundtruth to compare with the moduli values derived
from the FWD data and the ones suggested by the laboratory tests.
INPUTS TO THE PROGRAMS
General Inputs for All 4 Programs
Thickness of Layers
All of the pavement structures under study have been made up of an
asphalt surface, on a crushed stone base. These two layers rest on a subgrade
whose type is defined by mechanical tests. The thickness of each layer
at each station has also been measured.
Moduli Range and Poisson’s Ratio
In the backcalculation process (Table 1) from SHRP backcalculation
report was used to define the surface and base layers moduli ranges. For the
subgrade, the measured CBR values were used to estimate the resilient modulus
using the relationships given in Table 2 and Fig.
4 from Iran Highway Asphalt Paving Code (2000).
Based on the assembled data the following ranges were thought reasonable
for the modulus of subgrade for each site:
• 
Zanjan Tabriz: 70 to 105 MPa 
• 
Eivanekey Garmsar and Garmsra Semnan: 70 to 140 MPa 
• 
Rafsanjan Airport: 170 to 210 MPa 
Additional Inputs to DBSID
Unit Weight (γ)
The value of 2250 kg m^{3} was selected for all asphalt concrete
layers. For base and subgrade Table 3 was used for determining
γ, if it was not clearly mentioned in the soil mechanic test results. Lastly,
the following values for the unit weight were selected:
• 
ZanjanTabriz: 1600 kg m^{3} for the subgrade, 1902 kg
m^{3} for the base layer 
• 
EivanekeyGarmsar and GarmsraSemnan: 1920 kg m^{3} for
both base and subgrade layers 
• 
Rafsanjan Airport: 1920 kg m^{3} for both base and subgrade
layers 
Table 1: 
Default moduli ranges and Poisson’s Ratio values (William,
1999) 

Table 2: 
Relationship between CBR and resilient modulus for the CBR
values less than 25 


Fig. 4: 
Relationship between CBR and resilient modulus for the CBR
values greater than 25 
Table 3: 
Table of typical m values and unit weight for base course
and subgrade materials (Lytton, 1989) 

Damping Ratio
From the sensitivity analysis performed it was found that the damping ratio
does not significantly effect the deflection time history a lot within the range
of this parameter (Fernando and Liu, 1993). The suggested
default value of 4% was fixed in the DBSID program.
After the completion of the input data, the programs were run and the
outputs were reviewed for reasonableness.
DISCUSSION
The elastostatic programs terminate when the convergence condition between
deflections is satisfied. The Root Mean Squared Errors (RMSE’s)
higher than 5% were not accepted. Reviewing the results is always necessary
because every predicted modulus may not be acceptable. For example, the
moduli which are equal to the lower or upper limit of the acceptable range
are difficult to justify given the known quality of the materials.
In EVERCALC program, it was found that moduli convergence usually happened
earlier than deflection convergence. So, the program was rerun by changing
the moduli range to achieve the best deflection fit.
The DBSID program tries to create the best fit between measured and computed
deflection time histories. Figure 5a and b present the results from this
study. The output files contain the RMSE between all deflection values
(during the 60 msec load duration) and also the errors between the peak
values of calculated and measured deflections at each sensor. The surface
moduli are most dependent on the deflection measured at the sensor in
the middle of the load plate (W1 at 0 inches), so in the final evaluation,
not only the RMSE have been considered, but also the peak error at the
W1 sensor.


Fig. 5: 
Measured and calculated deflections at one sample station
(output of DBSID) (a) Low RMS and peak errors and (b) Low RMS but high peak
errors 
In most of the DBSID output files from the different sections on the
ZanjanTabriz highway both RMS and peak errors were high. The time histories
of the measured deflections also had irregular shapes. The cause of this
was not determined in this study but it could be related to the following
reasons:
• 
Equipment Problems: FWD measurements were not taken with good accuracy 
• 
The pavement was distressed at the time of testing and the surface
and base condition may have impacted sensor readings 
These output files were excluded from the final conclusion.
CLASSIFICATION OF RESULTS
After the termination of every run and omitting unacceptable values of
moduli (the moduli which are not in the defined range) for each deflection
input file, the results were classified in two ways:
• 
Scatter charts such as the one shown in Fig. 6
were developed to monitor the variation of modulus in consecutive
stations 
• 
Bar charts such as the one shown in Fig. 7 compare
the average moduli of different layers 

Fig. 6: 
Moduli changes in different stations. This example result
is due to Apron number 3 of Rafsanjan airport (a) Surface (b) Base and (c)
Subgard 

Fig. 7: 
Comparison of the average moduli. This example result is due
to Apron number 3 of Rafsanjan airport (a) Surface and (b) Base and Subgrade 
Note that the average values in Fig. 7 are obtained
from scattered values of Fig. 6.
CONCLUSIONS
FWD data collected from four different test sites were analyzed by three
elastostatic backcalculation programs: MODULUS6.0, ELMOD5.0 and EVERCALC5.0
and a dynamic backcalculation program: DBSID. Qualitative comparison of
the predicted moduli are presented in the following two sections:
Comparison of Three ElastoStatic Programs
• 
The moduli change in consecutive stations and the average moduli
of surface layer, from the outputs of all the three programs show
good consistency in all three cases 
• 
ELMOD5.0 overestimates the resilient moduli of the subgrade in some
cases 
• 
MODULUS6.0 and EVERCALC5.0 show good consistency in results for
all three layers 
• 
As the depth to bedrock is changeable in MODULUS6.0, if the reasonable
range of the subgrade moduli is clear, the user can alter the depth
to bedrock until the desired values of subgrade moduli are obtained.
This is also possible in EVERCALC5.0, but the process is faster in
MODULUS6.0 
• 
Coefficients of Variance (CV), which is computed by dividing the
variance by the mean value, for the predicted subgrade moduli of MODULUS6.0
are lower than that of the other programs (not greater than 19%) 
• 
As the backcalculation process often starts from predicting the
subgrade modulus, a proper estimate of this value is critical for
the rest of the process. The subgrade moduli of MODULUS 6.0 were in
the range anticipated from known site conditions and from the CBR
data described earlier 
Comparison of the ElastoStatic Analysis of MODULUS6.0 with the Dynamic
Analysis of DBSID
• 
Subgrade moduli predicted by DBSID are almost 2.5 times greater
than subgrade moduli predicted by MODULUS 6.0. The ratio of the depth
of stiff layer from DBSID to MODULUS 6.0 is almost the same value
(2.5) 
• 
Base moduli calculated by MODULUS6.0 are almost 1.3 times greater
than that of DBSID. So, the base moduli of DBSID are changeable to
the base moduli of DBSID by a factor of 0.77 
• 
DBSID’s surface moduli are often equal or greater (up to 20%)
than that of MODULUS 6.0 
• 
Dynamic backcalculation is a time consuming task, analyzing each
station by DBSID takes about 10 min by normal PCS 
So, by considering all the mentioned points, MODULUS6.0 is suggested
as the most appropriate program for the cases under study.
ACKNOWLEDGMENT
The authors are grateful to Dr. Emanuel Fernando of Texas A and M University
for providing the DBSID program.