BACKCALCULATION OF ASPHALT
CONCRETE COMPLEX MODULUS CURVE
BY LAYERED VISCOELASTIC SOLUTION
By
Ligang Lei

A DISSERTATION
Submitted to
Michigan State University
In partial fulfillment of the requirements
For the degree of
DOCTOR OF PHILOSOPHY

Civil Engineering
2011

ABSTRACT
BACKCALCULATION OF ASPHALT CONCRETE COMPLEX MODULUS CURVE BY
LAYERED VISCOELASTIC SOLUTION
By
Ligang Lei
In-situ evaluation of material properties is very important for estimating the structural
adequacy of a pavement section under traffic loading. The fundamental material property for the
AC layer in a flexible pavement is the complex modulus |E*| or the equivalent relaxation
modulus E(t). In-situ |E*| can be used as a quality control tool as well as in estimation of
remaining service life of existing pavements, and it is a critical input to the new
Mechanism-Empirical Pavement Design Guide (M-E PDG). In this research, a method is
presented to backcalculate |E*| or E(t) for the AC layer, in addition to base and subgrade layer
moduli, using Falling Weight Deflectometer (FWD) surface deflection time histories. First, the
deflection time histories from FWD test data were separated as the dynamic response part, and
the viscoelastic response part: The time delay of each sensor is related to wave propagation
through the pavement system while the shifted deflection time history without time delay can be
considered as the viscoelastic response of the pavement. Second, the time delay of each sensor
was used to estimate the elastic modulus of the subgrade layer, based on the wave propagation
theory. Third, a forward layered viscoelastic solution was developed based on Schapery’s
‘quasi-elastic’ approximation. Finally, backcalculation was numerically done by Newton’s
method, and followed by the MATLAB internal function ‘fminsearch’, to match the predicted
deflection time histories from estimated modulus of each layer against the shifted measurement
from the FWD test. The backcalculation results show very good agreement with the actual

modulus values for both numerical examples and field FWD test data. The error in the base and
subgrade moduli is generally less than 3%. The AC layer relaxation and complex modulus curves
match the actual functions within the ranges of 0.0001 to 0.1 sec, and 10 to 10,000 Hz,
respectively. The sensitivity analysis shows that errors in the deflection time-histories are the
most significant factor affecting the backcalculated AC relaxation and complex modulus curves.
Also, the layer thicknesses should be as accurate as possible.

TO MY PARENTS, FAMILY AND FRIENDS

iv

ACKNOWLEDGEMENT
The author would like to express his great appreciation to his major advisor Dr. Karim Chatti,
Professor of Civil Engineering at Michigan State University, for his continuous academic and
financial support throughout the author’s study at Michigan State University and while preparing
this dissertation. The author would like to thank his co-advisor, Dr. Emin Kutay, Assistant
Professor of Civil Engineering, for his valuable suggestions and kind help. This dissertation
could not have been completed without their help. The author also appreciated help from other
committee members, Dr. Gilbert Baladi, Professor of Civil Engineering, and Dr. Tom Pence,
Professor of Mechanical Engineering, for their continuous advice and support.

Thanks are also due to the Federal Highway Administration (FHWA) for financially sponsoring
this research, and to the author’s graduate student colleagues at the Department of Civil and
Environmental Engineering at MSU.

Last but not least, the author would like to thank his wife, Huajing Yang, for her continuous and
unselfish support and company throughout his stay at Michigan State University.

v

TABLE OF CONTENTS

LIST OF TABLES......................................................................................................................... ix
LIST OF FIGURES..................................................................................................................... xii
Chapter 1. Introduction ................................................................................................................. 1
Chapter 2. Literature Review and Background ............................................................................ 5
2.1 Backcalculation Approach ................................................................................................. 5
2.1.1 Static Backcalculation Methods..................................................................................... 6
2.1.2 Dynamic Backcalculation Methods ............................................................................... 7
2.2 Primary Response Behavior of Asphalt Concrete............................................................ 8
2.3 Interconversion between E(t) and |E*| for AC Mixture ................................................ 13
2.4 Interconversion between E(t) and D(t) for AC Mixture ................................................ 15
2.5 Motivation for the Study .................................................................................................. 19
Chapter 3. Estimating the Subgrade Elastic Modulus using Wave Propagation Method ........ 22
3.1 Introduction....................................................................................................................... 22
3.2 Proposed Method to Estimate k-value ............................................................................ 23
3.2.1 Theoretical Formulation for Spectral Element Method ............................................... 24
3.2.2 Dominant Frequency of the FWD Test ........................................................................ 31
3.2.3 Numerical Example to Calculate k-value .................................................................... 32
3.2.4 Sensitivity Analysis for k-value ................................................................................... 33
3.3 Proposed Procedure to Estimate Subgrade Elastic Modulus ....................................... 38
3.4 Verification of the Proposed Method............................................................................... 42
3.4.1 SAPSI........................................................................................................................... 42
3.4.2 LAMDA....................................................................................................................... 43
3.4.3 Case of Subgrade without Ground Water Table........................................................... 44
3.4.4 Case of Subgrade with GWT ....................................................................................... 46
3.5 Validating the Proposed Method Using Case Studies from LTPP................................ 49
3.5.1 Example of FWD Test at LTPP 04-1036 Station in 1998 ............................................ 49
3.5.2 Example of FWD Test at LTPP 32-0101 Station in 1996 ............................................ 52
3.5.3 Example of the FWD Test at LTPP Station 06-0565 in 1999 ...................................... 56
3.5.4 Discussion of Case Studies Obtained from LTPP........................................................ 60
3.5.4.1 Effect of Confinement....................................................................................... 60
vi

3.5.4.2 Effect of Sampling Intervals ............................................................................. 63
3.6 Numerical Example of Bedrock under Subgrade layer ................................................ 64
3.7 Summary............................................................................................................................ 67
Chapter 4. Layered Viscoelastic Forward Solution .................................................................... 68
4.1 Introduction....................................................................................................................... 68
4.2 Layered Viscoelastic Forward Solution Algorithm........................................................ 68
4.2.1 Theoretical Background............................................................................................... 68
4.2.2 Algorithm Steps ........................................................................................................... 71
4.3 Verification of the Layered Viscoelastic Forward Solution........................................... 76
4.3.1 Verification of Quasi-Elastic Approximation............................................................... 76
4.3.2 Comparison between the Layered VE Solution and Dynamic Solutions .................... 78
4.3.3 Verification Examples using Pavement Structures from the LTPP Database .............. 83
4.3.4 A Numerical Example with the Bedrock ..................................................................... 89
4.4 Time Efficiency for Layered VE Forward Solution....................................................... 91
4.5 Summary............................................................................................................................ 93
Chapter 5. Backcalculation Algorithm........................................................................................ 94
5.1 Introduction....................................................................................................................... 94
5.2 Variable Identification ...................................................................................................... 94
5.3 Variable Boundaries.......................................................................................................... 95
5.3.1 Boundary Limits for |E*|.............................................................................................. 95
5.3.2 Boundary Limits for E(t) ............................................................................................. 97
5.4 Inverse Solution Theory ................................................................................................... 99
5.4.1 Basic Theory Background............................................................................................ 99
5.4.2 SVD Method .............................................................................................................. 102
5.4.3 Truncating Singular Values ........................................................................................ 103
5.5 VE backcalculation of E(t) ............................................................................................. 104
5.6 VE backcalculation of |E*| ............................................................................................. 107
5.6.1 Forward calculation of E(t) from |E*| and φ .............................................................. 107
5.6.2 Calculation of E(t) and φ from |E*| by Iteration ........................................................ 108
5.6.3 Calculation of E(t) from |E*(ω)| by Iteration φ(ω) .....................................................112
5.6.4 Flowchart of VE backcalculation of |E*| ....................................................................114
5.7 One Numerical Example of Backcalculation Result.....................................................116
5.7.1 Result of VE backcalculation of E(t) ..........................................................................116
5.7.2 Result of VE backcalculation of |E*| ..........................................................................118
5.8 Summary.......................................................................................................................... 120
Chapter 6. Verification and Sensitivity Analysis ....................................................................... 121
vii

6.1 Introduction..................................................................................................................... 121
6.2 Verification of the VE Backcalculation Algorithm ...................................................... 121
6.2.1 VE Backcalculation using Time Histories from VE Solution.................................... 123
6.2.1.1 Numerical Example 1 ..................................................................................... 124
6.2.1.2 Numerical Example 2 ..................................................................................... 129
6.2.1.3 Numerical Example 3 ..................................................................................... 133
6.2.1.4 Numerical Example 4 ..................................................................................... 137
6.2.1.5 Discussion ....................................................................................................... 141
6.2.2 VE Backcalculation using Time Histories from Dynamic Solutions......................... 141
6.2.2.1 Example 5: Michigan Site............................................................................... 141
6.2.2.2 Example 6: Texas Site..................................................................................... 145
6.2.2.3 Example 7: Florence Site ................................................................................ 150
6.2.2.4 Example 8: Kansas Site .................................................................................. 154
6.2.2.5 Discussion ....................................................................................................... 158
6.2.3 VE Backcalculation using Field FWD Test Data....................................................... 158
6.2.3.1 FWD Test Conducted at L9S4 ........................................................................ 159
6.2.3.2 FWD Test Conducted at L10S4 ...................................................................... 161
6.2.3.3 FWD Test Conducted at L11S4 ...................................................................... 163
6.2.3.4 Discussion ....................................................................................................... 165
6.3 Sensitivity Analysis for the AC Layer ........................................................................... 165
6.3.1 Sensitivity Analysis for Deflection Measurement of the First Sensor....................... 168
6.3.2 Sensitivity Analysis for Time Shift Error of the First Sensor .................................... 172
6.3.3 Sensitivity Analysis for Post-Peak Time Offset of the First Sensor .......................... 177
6.3.4 Sensitivity Analysis for Layer Thickness .................................................................. 181
6.3.5 Discussion of Sensitivity Analyses ............................................................................ 186
6.4 Summary.......................................................................................................................... 186
Chapter 7. Conclusions and Recommendations ....................................................................... 188
7.1 Summary.......................................................................................................................... 188
7.2 Conclusions...................................................................................................................... 189
7.3 Recommendations ........................................................................................................... 190
BIBLIOGRAPHY....................................................................................................................... 192

viii

LIST OF TABLES

Table 1. Prony series coefficients for E(t) and D(t) for a typical AC mixture ...................... 18
Table 2: The equivalent frequency for different load duration times.................................... 32
Table 3: The average k-value for different cases of Poisson’s ratio ..................................... 37
Table 4: Basic information of a three-layer pavement structure without GWT ................... 44
Table 5: Time delay of each sensor simulated by SAPSI ..................................................... 44
Table 6: Basic information for a three-layer pavement structure with GWT ....................... 46
Table 7: Time delay of each sensor by LAMDA .................................................................. 47
Table 8: Time delay and peak deflection of each sensor due to the FWD test at station
04-1036 in 1998 ............................................................................................................ 50
Table 9: Structural information and MODCOMP5 backcalculated modulus at LTPP station
04-1036 in 1998 ............................................................................................................ 50
Table 10: Time delay and maximum deflection of sensors during the FWD test at LTPP
station 32-0101 in 1996 ................................................................................................ 53
Table 11: Basic information and MODCOMP5 backcalculated modulus at LTPP station
32-0101 in 1996 (Trial: I) ............................................................................................. 54
Table 12: Basic information and MODCOMP5 backcalculated modulus at LTPP station
32-0101 in 1996 (Trial: II) ............................................................................................ 54
Table 13: Elastic modulus of the subgrade layer at LTPP station 32-0101 in 1996 ............. 55
Table 14: Time delay and maximum deflection of the sensors due to the FWD test at LTPP
station 06-0565 in 1999 ................................................................................................ 57
Table 15: Structural information and MODCOMP5 backcalculated modulus at LTPP station
06-0565 in 1999 (Trial: I) ............................................................................................. 58
Table 16: Structural information and MODCOMP5 backcalculated modulus at LTPP station
06-0565 in 1999 (Trial: II) ............................................................................................ 58
Table 17: Elastic modulus of the subgrade layer at station 06-0565 in 1999 ....................... 59
Table 18: Results of data regression of resilient modulus on stress parameters using Uzan
ix

and Cornell Models....................................................................................................... 61
Table 19: The stress status of the subgrade layer in the FWD test for different K0 values... 62
Table 20: The range of elastic moduli predicted by nonlinear models, measured in the lab,
estimated by the proposed method, and backcalculated by MODCOMP5 for the case
studies (MPa) ................................................................................................................ 63
Table 21: Arrival time for S-wave and R-wave for each sensor due to the FWD test.......... 64
Table 22: Material properties of a three-layer pavement structure with bedrock ................. 65
Table 23: Time delay of each sensor by SPASI simulation with bedrock ............................ 66
Table 24: Material properties of the pavement structure used in the verification example .. 76
Table 25: Properties of a pavement used in comparison of dynamic and VE solutions ....... 80
Table 26: Several examples of pavement structure in the SPS-1 project ............................. 83
Table 27: Layer information for each pavement structure.................................................... 84
Table 28: The computation times for different numbers of discrete time steps.................... 92
Table 29: The |E*| parameters for different AC mixtures ..................................................... 96
Table 30: The average and boundaries for the four parameters of |E*|................................. 97
Table 31: The E(t) parameters for different AC mixtures ..................................................... 98
Table 32: The mean and boundaries for the four parameters of E(t) .................................... 99
Table 33: Basic information of a three-layer pavement structure and backcalculated result
......................................................................................................................................116
Table 34：Basic information of pavement structure and backcalculated result for numerical
example 1 .................................................................................................................... 124
Table 35：Basic information of pavement structure and backcalculated result for numerical
example 2 .................................................................................................................... 129
Table 36：Basic information of pavement structure and backcalculated result for numerical
example 3 .................................................................................................................... 133
Table 37：Basic information of pavement structure and backcalculated result for numerical
example 4 .................................................................................................................... 137
Table 38：Basic information of Michigan pavement structure and backcalculated results142
x

Table 39：Basic information of Texas pavement structure and backcalculated results..... 146
Table 40: Basic information of a Florence pavement structure and backcalculated result. 150
Table 41: Basic information of a Kansas pavement structure and backcalculated result ... 154
Table 42: Basic information of L9S4 pavement structure and backcalculated result......... 159
Table 43: Basic information of L10S4 pavement structure and backcalculated result....... 161
Table 44: Basic information of L11S4 pavement structure and backcalculated result ....... 163
Table 45: Layer information of each pavement structure for sensitivity analysis .............. 166

xi

LIST OF FIGURES

Figure 1: Organization of this dissertation.............................................................................. 3
Figure 2: Sinusoidal Stress and Strain in cyclic loading......................................................... 9
Figure 3: Mastercurve |E*| of one asphalt concrete mixture (PG64-28) .............................. 10
Figure 4: Phase angle φ(ω) of one asphalt concrete mixture (PG64-28)...............................11
Figure 5: Physical meaning of each parameter on the sigmoidal function ........................... 12
Figure 6: Relaxation modulus |E(t)| of one asphalt concrete mixture (PG64-28)................. 13
Figure 7: The generalized Maxwell model ........................................................................... 14
Figure 8: E(t) and D(t) from interconversion for a typical AC mixture................................ 19
Figure 9: Wave system from surface point source in ideal medium ..................................... 30
Figure 10: Phase angle of each point at the surface for a specific example (t =0 s)............. 32
Figure 11: Sensitivity analysis for k-value (ν= 0.20)............................................................ 34
Figure 12: Sensitivity analysis for k-value (ν= 0.25)............................................................ 34
Figure 13: Sensitivity analysis for k-value (ν= 0.30)............................................................ 35
Figure 14: Sensitivity analysis for k-value (ν= 0.35)............................................................ 35
Figure 15: Sensitivity analysis for k-value (ν= 0.40)............................................................ 36
Figure 16: Sensitivity analysis for k-value (ν= 0.45)............................................................ 36
Figure 17: Sensitivity analysis for k-value (ν= 0.49)............................................................ 37
Figure 18: Example of FWD load and deflection time histories at LTPP station 04-1036 in
1998............................................................................................................................... 39
Figure 19: Time delay vs. location for the example of FWD test at LTPP station 04-1036 in
1998............................................................................................................................... 40
Figure 20: Flow chart to estimate subgrade elastic modulus using k-value ......................... 41
Figure 21: Time delay vs. location for a numerical example of FWD test (without GWT). 45
Figure 22: LAMDA simulated time-history due to the FWD test (with GWT) ................... 47
xii

Figure 23: LAMDA simulated time-delay for the FWD test (with GWT) ........................... 48
Figure 24: Deflection basin from FWD test and MODCOMP5 Simulation at LTPP station
04-1036 in 1998 ............................................................................................................ 51
Figure 25: Example of FWD test results at LTPP station 32-0101 in 1996.......................... 53
Figure 26: Deflection basin due to the FWD test and two cases of MODCOMP5 simulation
at LTPP station 32-0101 in 1996................................................................................... 55
Figure 27: Example of the FWD test history at LTPP station 06-0565 in 1999 ................... 57
Figure 28: Deflection basin during the FWD test and two cases of MODCOMP5 simulation
at LTPP station 06-0565 in 1999................................................................................... 59
Figure 29: SAPSI simulated time-history of FWD test with bedrock underneath ............... 65
Figure 30: Time delay vs. location for a numerical example of FWD test with bedrock..... 66
Figure 31: Typical geometry of a pavement structure .......................................................... 70
Figure 32: Discretization of stress history ............................................................................ 72
Figure 33: Discretization of the relaxation modulus mastercurve ........................................ 72
Figure 34: Deflections calculated for points at different distances from the centerline of the
circular load at the surface ............................................................................................ 73
Figure 35: Illustration of the

dσ (τ j )

in Eq. (32b) for each time step

τ j .................... 74

Figure 36: Example of a pavement structure used for VE forward calculation.................... 75
Figure 37: Examples of computed VE surface deflections at different radial distances from
the centerline of the load............................................................................................... 75
Figure 38: Implemented and actual E(t) for one AC mixture ............................................... 77
Figure 39: Implemented and actual E(t) for one AC mixture ............................................... 77
Figure 40: FWD test simulation by layered VE solution for a half-space AC layer ............ 79
Figure 41: Sensor deflection time histories predicted by SAPSI.......................................... 80
Figure 42: Sensor deflection time histories predicted by LAMDA...................................... 81
Figure 43: Comparison of deflection time histories from VE and SAPSI solutions ............ 82
Figure 44: Comparison of deflection time histories from VE and LAMDA solutions......... 82
xiii

Figure 45: Time-shifted dynamic solution and VE solution for case 116............................. 85
Figure 46: Time-shifted dynamic solution and VE solution for case 117............................. 86
Figure 47: Time-shifted dynamic solution and VE solution for case 120 ............................ 87
Figure 48: Time-shifted dynamic solution and VE solution for case 123 ............................ 88
Figure 49: Time-shifted SAPSI and VE solution for shallow bedrock................................. 90
Figure 50: Sensitivity analysis for relative error of maximum deflection between VE
solution and SAPSI with bedrock depth ....................................................................... 91
Figure 51: Random Function to avoid local minimum problem......................................... 105
Figure 52: Algorithm of backcalculation procedure of E(t)................................................ 106
Figure 53: Algorithm for Calculation of E(t) and φ from |E*| by iteration..........................111
Figure 54: Algorithm for Calculation of E(t) from |E*| by iteration φ ................................113
Figure 55: |E*| input and calculated after converting to E(t) for one AC mixture...............114
Figure 56: Algorithm of backcalculation procedure of |E*|.................................................115
Figure 57: Input and backcalculated time-histories of E(t) in the numerical FWD test ......117
Figure 58: Input and backcalculated E(t) for the numerical FWD test................................117
Figure 59: Input and backcalculated time-histories of |E*| in the numerical FWD test ......119
Figure 60: Input and backcalculated |E*| for the numerical FWD test ................................119
Figure 61: Available range in E(t) function from field FWD test backcalculation............. 122
Figure 62: Useful range in |E*| function from field FWD test backcalculation ................. 123
Figure 63: Input and simulated time-histories by backcalculation of E(t) in example 1.... 125
Figure 64：Input and backcalculated E(t) in example 1 .................................................... 125
Figure 65: Input and simulated time-histories by backcalculation of |E*| in example 1 .... 126
Figure 66: Input and backcalculated |E*| in example 1 ...................................................... 127
Figure 67: Input and backcalculated φ in example 1 .......................................................... 128
Figure 68: Input and simulated time-histories by backcalculation of E(t) in example 2.... 130

xiv

Figure 69：Input and backcalculated of E(t) in example 2................................................ 131
Figure 70: Input and simulated time-histories by backcalculation of |E*| in example 2 .... 132
Figure 71: Input and backcalculated of |E*| in example 2 .................................................. 132
Figure 72: Input and simulated time-histories by backcalculation of E(t) in example 3.... 134
Figure 73：Input and backcalculated E(t) in example 3 .................................................... 135
Figure 74: Input and simulated time-histories by backcalculation of |E*| in example 3 .... 135
Figure 75: Input and backcalculated |E*| in example 3 ...................................................... 136
Figure 76: Input and simulated time-histories by backcalculation of E(t) in example 4.... 138
Figure 77：Input and backcalculated E(t) in example 4 .................................................... 138
Figure 78: Input and simulated time-histories by backcalculation of |E*| in example 4 .... 139
Figure 79: Input and backcalculated |E*| in example 4 ...................................................... 140
Figure 80: Input and simulated time-histories by backcalculation of E(t) in Michigan site143
Figure 81：Input and backcalculated E(t) in Michigan site............................................... 143
Figure 82: Input and simulated time-histories by backcalculation of |E*| in Michigan site144
Figure 83: Input and backcalculated |E*| in Michigan site ................................................. 145
Figure 84：Input and simulated time-histories by backcalculation of E(t) in Texas site... 147
Figure 85：Input and backcalculated E(t) in Texas site ..................................................... 147
Figure 86：Input and simulated time-histories by backcalculation of |E*| in Texas site ... 148
Figure 87：Input and backcalculated |E*| in Texas site ..................................................... 149
Figure 88: Input and backcalculation simulated time-histories using E(t) in Florence site 151
Figure 89: Input and backcalculated E(t) in Florence site .................................................. 151
Figure 90: Input and backcalculation simulated time-histories using |E*| in Florence site 152
Figure 91: Input and backcalculated |E*| in Florence site .................................................. 153
Figure 92: Input and simulated time-histories by backcalculation of E(t) in Kansas site .. 155
Figure 93: Input and backcalculated E(t) in Kansas site..................................................... 155
xv

Figure 94: Input and simulated time-histories by backcalculation of |E*| in Kansas site .. 156
Figure 95: Input and backcalculated |E*| in Kansas site..................................................... 157
Figure 96: Input and simulated time-histories by backcalculation of E(t) in L9S4 site ..... 160
Figure 97: Input and backcalculated E(t) in L9S4 site ....................................................... 160
Figure 98: Input and simulated time-histories by backcalculation of E(t) in L10S4 site ... 162
Figure 99: Input and backcalculated E(t) in L10S4 site ..................................................... 162
Figure 100: Input and simulated time-histories by backcalculation of E(t) in L11S4 site . 164
Figure 101: Input and backcalculated E(t) in L11S4 site.................................................... 164
Figure 102: Error caused by time shift of entire pulse for one sensor................................ 167
Figure 103: Error caused by expansion/contraction of deflection pulse in the post-peak zone
for one sensor.............................................................................................................. 167
Figure 104: Sensitivity analysis of |E*| for deflection measurement error (Case A).......... 168
Figure 105: Sensitivity analysis of E(t) for deflection measurement error (Case A).......... 169
Figure 106: Sensitivity analysis of |E*| for deflection measurement error (Case B).......... 169
Figure 107: Sensitivity analysis of E(t) for deflection measurement error (Case B).......... 170
Figure 108: Sensitivity analysis of |E*| for deflection measurement error (Case C).......... 170
Figure 109: Sensitivity analysis of E(t) for deflection measurement error (Case C).......... 171
Figure 110: Sensitivity analysis of |E*| for deflection measurement error (Case D).......... 171
Figure 111: Sensitivity analysis of E(t) for deflection measurement error (Case D).......... 172
Figure 112: Sensitivity analysis of |E*| for time shift error of sensor 1 (Case A) .............. 173
Figure 113: Sensitivity analysis of E(t) for time shift error of sensor 1 (Case A) .............. 173
Figure 114: Sensitivity analysis of |E*| for time shift error of sensor 1 (Case B) .............. 174
Figure 115: Sensitivity analysis of E(t) for time shift error of sensor 1 (Case B) .............. 174
Figure 116: Sensitivity analysis of |E*| for time shift error of sensor 1 (Case C) .............. 175
Figure 117: Sensitivity analysis of E(t) for time shift error of sensor 1 (Case C) .............. 175
Figure 118: Sensitivity analysis of |E*| for time shift error of sensor 1 (Case D) .............. 176
xvi

Figure 119: Sensitivity analysis of E(t) for time shift error of sensor 1 (Case D) .............. 176
Figure 120: Sensitivity analysis of |E*| for post-peak time offset of sensor 1 (Case A)..... 177
Figure 121: Sensitivity analysis of E(t) for post-peak time offset of sensor 1 (Case A) .... 178
Figure 122: Sensitivity analysis of |E*| for post-peak time offset of sensor 1 (Case B)..... 178
Figure 123: Sensitivity analysis of E(t) for post-peak time offset of sensor 1 (Case B) .... 179
Figure 124: Sensitivity analysis of |E*| for post-peak time offset of sensor 1 (Case C)..... 179
Figure 125: Sensitivity analysis of (t) for post-peak time offset of sensor 1 (Case C)....... 180
Figure 126: Sensitivity analysis of |E*| for post-peak time offset of sensor 1 (Case D) .... 180
Figure 127: Sensitivity analysis of E(t) for post-peak time offset of sensor 1 (Case D) .... 181
Figure 128: Sensitivity analysis of |E*| for layer thickness (Case A) ................................. 182
Figure 129: Sensitivity analysis of E(t) for layer thickness (Case A)................................. 182
Figure 130: Sensitivity analysis of |E*| for layer thickness (Case B) ................................. 183
Figure 131: Sensitivity analysis of E(t) for layer thickness (Case B)................................. 183
Figure 132: Sensitivity analysis of |E*| for layer thickness (Case C) ................................. 184
Figure 133: Sensitivity analysis of E(t) for layer thickness (Case C)................................. 184
Figure 134: Sensitivity analysis of |E*| for layer thickness (Case D)................................. 185
Figure 135: Sensitivity analysis of E(t) for layer thickness (Case D)................................. 185

xvii

Chapter 1. Introduction
Asphalt concrete (AC) mixtures consist of asphalt binder and aggregates from crushed stone
or from natural resources. They were first used for pavement in 1870 in Newark, New Jersey.
According to the Federal Highway Administration (FHWA, 2001), 2.5 million miles of road are
paved in the USA, and 94% are paved with AC, i.e., about 2.3 million miles of road are
surface-covered by AC, or approximately 9 times the distance to the moon. About 300 million
tons of AC mixtures are used in highway construction every year (Papagiannakis and Masad,
2008).
Before the 1920s, only the thickness of pavement was considered in pavement design and
analysis. It was believed the thicker the pavement, the longer the pavement would last. After
many years of experience, three key factors are found to significantly affect pavement life: traffic
or loading, environment and material properties. Traffic information is estimated by regional
planning departments, while national weather stations can provide environmental conditions for
a specific location. The remaining challenge for pavement engineering is the material properties
of the pavement, especially the onsite materials used many years ago.
In recent years, the falling weight deflectometer (FWD) test has been the typical method for
evaluating the material properties of in-service pavement. In order to accurately evaluate
pavement structure, several FWD tests with different load levels are usually performed at the
same location. The structural condition, or layer modulus, is the key factor for determining
pavement rehabilitation strategies.
The material property is usually obtained through static or dynamic backcalculation methods,
1

in which layer moduli are determined by matching the deflection basin measured under a known
load in the FWD test, with the deflection basin generated through a theoretical model of the
pavement (Ji, 2005).
Static backcalculation is mainly based on the layered linear elastic solution, while dynamic
backcalculation modifies the linear elastic property into the damped elastic property for each
physical layer. However, both models of the pavement structure are theoretically incorrect,
because the primary response, i.e., small load-level response, of AC mixture is viscoelastic,
rather than elastic. The mechanical behaviors of viscoelastic materials are quite different from
those of the elastic or damped elastic material, for example, granular materials, aggregate base,
or subgrade.
The objective of this research is to develop a robust methodology to backcalculate the |E*|
and E(t) mastercurve of the AC layer, and the elastic moduli of unbound layers.
The dynamic effect due to the FWD test should be investigated first, and the pavement
response is expected to be simplified as a viscoelastic response. Additionally, the elastic modulus
of the subgrade layer can be estimated or identified from the dynamic response. After that, a
pavement structure is simulated as viscoelastic system, and backcalculation is done in
viscoelastic solution. A brief flowchart and the explanation of this dissertation are shown as
following:

2

Chapter 3: Wave propagation
and estimation Esubgrade
AC
AC

Ε(t), ν1
Ε ,ν

Base

H=Inf

2

Subgrade

Chapter 4: Layered viscoelastic
forward solution

2

Chapter 5: Backcalculation
of E(t)/ E* and Eunbound

Εn, νn

Chapter 6: Backcalculation
results & sensitivity analysis

Figure 1: Organization of this dissertation
(For interpretation of the references to color in this and all other figures, the reader is
referred to the electronic version of this dissertation.)
Chapter 2 provides a review of relevant literature on the backcalculation approaches from
field FWD tests, background on the material properties of the AC mixture as viscoelastic
materials, and the interconversions between |E*|, E(t) and D(t) within the frame work of
viscoelasticity. The purpose of the study is briefly given at the end of this chapter.
Chapter 3 presents a method for calculation of Esubgrade based on wave propagation theory.
The concept of equivalent frequency is introduced, and the application of Rayleigh wave in
analysis of FWD test data is explained. As a validation, four cases of field tests from Long
Term Pavement Performance (LTPP) database are demonstrated. The merits and limitations
are given.

3

Chapter 4 introduces a layered viscoelastic forward solution. The response of viscoelastic
multi-layered pavement system under FWD test is presented. The pavement is modeled as a
system of the horizontal layers, and AC mixture is simulated as viscoelastic material, and the
rest layers are isotropic and linear elastic. The forward solution is verified against the
semi-analytical viscoelastic solution of a step surface load, and then checked against
well-known dynamic solution SAPSI and LAMDA.
Chapter 5 contains the details of the backcalculation algorithm, and some brief information
on the backcalculation algorithm developed in MATLAB. A hypothetical FWD test
deflection time histories, simulated by dynamic solution SAPSI, was used as the input
parameters, and the modulus of each physical layer was backcalculated to verify the
accuracy of the algorithm.
Chapter 6 presents the validation of the viscoelastic backcalculation algorithm using
numerical and field FWD test deflection time histories. Sensitivity analyses of various
parameters on the backcalculation were conducted. The effects of inaccuracies in deflections
and duration of deflection records were studied.
Chapter 7 contains a summary of the findings, the impact of the research, and
recommendations for future research.

4

Chapter 2. Literature Review and Background
2.1 Backcalculation Approach
The elastic moduli of pavement layers are obtained from backcalculation via field FWD tests.
To find the elastic modulus of each layer, backcalculation from FWD test data is usually carried
out by matching the measured deflections under a known load with theoretical deflections
generated by an analytical model of the pavement, usually by Newton’s method (Harichandran et
al, 1994). Such procedures usually use error minimization techniques to minimize either the
absolute or the square error, with or without weighting factors for sensors in the FWD test.
At present, pavement moduli can be backcalculated from the FWD deflection basin, using
the peak values of the deflection time-histories (static backcalculation) or using the FWD full
time-history (dynamic backcalculation). However, deflection basins under static loads differ
from those under dynamic or impulse loads because of geometric damping of pavement structure,
and dynamic effects, such as inertia, material damping, and resonance. Dynamic analysis, or the
time-history record, would therefore provide a more accurate estimation of the pavement moduli,
at the cost of time.
Additionally, the interpretation of data still remains problematic (Ji, 2005). This is due to the
limitation associated with the mechanical models incorporated into the backcalculation
procedures and the uniqueness of the inverse solutions. Nevertheless, during the past few
decades, there has been a significant improvement in the area of pavement modeling and
non-destructive test (NDT) techniques. In the following sub-sections, the static and dynamic

5

methods and backcalculation schemes are discussed.

2.1.1 Static Backcalculation Methods
The simplest way to model the behavior of flexible pavements is based on Boussinesq’s
(1885) solution that models a flexible pavement as a homogeneous, isotropic, and elastic
half-space. Later, Burmister (1943) modified this method into a two-layer system. Acum and Fox
(1951) advanced this model into a three-layer system. Newer methods, including CHEVRON
(Warren and Dieckmann, 1963) and KENLAYER (Huang, 2004) were implemented to analyze
the interface conditions between layers. Finite element methods, such as MICH-PAVE (Yeh,
1989) were also developed. The FWD test load is simulated as a static load with the maximum
magnitude of the impulse. The most recent layered elastic model is CHEVLAY2 (Stubstad et al,
2007).
Backcalculation programs typically perform on various forward computations to match the
deflection basins from FWD test to the computed deflections. There are three major groups of
static backcalculation methods, and each utilizes different techniques to reach the solution.
The first group is based on iteration techniques, where the layer moduli are repeatedly
adjusted until a suitable match between the calculated and measured deflection basins is obtained.
The typical programs for this method are MODCOMP (Irwin, 1994), BISDEF (Bush, 1985),
BOUSDEF (Roesset et al, 1995), and CHEVDEF (Bush and Alexander, 1985).
The second group is based on searching a database of deflection basins. A forward
calculation scheme is used to generate a database, which is then searched to find a best match for
the observed deflection basin. One typical example is MODULUS (Uzan, 1994).
The third group is based on the use of regression equations fitted to a database of deflection
basins generated by a forward calculation scheme. The LOADRATE program (Chua and Lytton,
6

1985) belongs to this category and uses regression equations generated from a database obtained
by using the ILLIPAVE (Raad and Figueroa, 1980) nonlinear finite element program.
A more detailed review on static backcalculation methods can be found elsewhere
(Mahmood, 1993).

2.1.2 Dynamic Backcalculation Methods
Most dynamic backcalculation methods use dynamic damped-elastic finite layer or finite
element models for their forward solutions. The finite layer solutions are based on Kausel’s
formulation (Kausel and Peek, 1982) which subdivides the medium into discrete layers that have
a linear displacement function in the vertical direction and satisfy the wave equation in the
horizontal direction. The solution is based on the premise that if the sub-layer thickness is small
relative to the wavelength of interest (around one tenth), it is possible to linearize the
transcendental functions and reduce them to algebraic expressions. The typical example of this
method is SAPSI (Chen, 1987). Al-Khoury et al (2001) developed an efficient forward solution,
LAMDA, for the dynamic analysis of flexible pavements using the spectral element technique
for the simulation of wave propagation in layered systems. The method is able to model each
layer as one element without the need for subdivision into several sub-layers. However, the
horizontal range (R) is required to simulate the vanishing of the wave toward the infinity.
Endiran (1999) developed a non-linear dynamic model in DYNARK, a computer program
accounting for non-linear properties in granular material as well as subgrade soil.
The backcalculation methods for dynamics are based on either frequency or time domain
solutions. In the frequency domain, the applied load and measured deflection time-histories are
transformed into the frequency domain by using the Fast Fourier Transform (FFT).
Backcalculation of layer parameters is done by matching the calculated steady-state (complex)
7

deflection basins with the frequency component of the measured sensor deflections at one or
more frequencies. In time domain backcalculation, the measured deflection time histories are
directly compared with the predicted results predicted by the forward program.
The advantages of backcalculation in the frequency domain are the computational efficiency
and theoretical validity. The disadvantage is that the results are very sensitive to truncation,
sampling time interval, and rest period after the loading impulse. Truncation of deflection history
of the sensors in the FWD tests is very common. The advantage of backcalculation in the time
domain is that the matching of deflection history can be achieved for any desired time interval.
The disadvantage is too many outputs in the forward calculation, so it is hard to converge for
forward and backcalculated results, perhaps even impossible. A more detailed review of dynamic
backcalculation can be found elsewhere (Ji, 2005).

2.2 Primary Response Behavior of Asphalt Concrete
It is well known that asphalt concrete (AC) behaves as a linear elastic material under small
strain conditions (Kim, 2008). A typical stress/strain behavior of AC in cyclic tests is shown in
Figure 2.

8

Stress/strain under cyclic loading

Normalized stress/strain

1

Normalized Stress
Normalized Strain

0.5

0

-0.5

-1
time

Figure 2: Sinusoidal Stress and Strain in cyclic loading
The time lag between stress and strain history illustrates that AC mixtures are viscoelastic
materials. Typically, the magnitude of dynamic complex modulus |E*| and the phase angle φ are
measured at different frequencies and at different loading temperatures. After eliminating the
effect of temperature by the time-temperature superposition principle, the mastercurve |E*| is
obtained, which is the magnitudes of dynamic modulus E* as a function of reduced frequency at
reference temperature. Experimental data show both |E*| and φ are only function of reduced
frequency ω for a specific mixture. One example of |E*| mastercurve is shown in Figure 3, and
the phase angle function is shown in Figure 4:

9

100000

|E*| (MPa)

10000

1000

E* Average

100

10
1.E-06

Sigmoid Fit, |E*| Mpa

1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Reduced Freq. (Hz)

Figure 3: Mastercurve |E*| of one asphalt concrete mixture (PG64-28)
(Kutay, 2008)

10

40.0

35.0

30.0

Phase Angle

25.0

20.0

y = 0.2186x3 - 0.53x2 - 6.1464x + 26.906
R2 = 0.9785

15.0

10.0

5.0

0.0
-6.0E+00

-5.0E+00

-4.0E+00

-3.0E+00

-2.0E+00

-1.0E+00

0.0E+00

1.0E+00

2.0E+00

3.0E+00

4.0E+00

Log fr

Figure 4: Phase angle φ(ω) of one asphalt concrete mixture (PG64-28)
(Kutay, 2008)
By regression, the mastercurve can be mathematically expressed as a sigmoidal function
with four parameters (a, b, c, d) obtained by data fitting.

log( E * ) = a +

b

1+ e

−c −d *log ( wR )

(1)

where: ωR is the reduced frequency at the reference temperature.
Each parameter in the sigmoidal function has a special physical meaning: ‘a’ represents
minimum modulus values at low reference frequency; ‘a + b’ indicates maximum modulus
values at high reference frequency; ‘c’ shows the horizontal position of the turning point; and ‘d’
influences the steepness of the function (rate of change between minimum and maximum
11

modulus values) (Kim, 2008). They are illustrated in Figure 5. The physical meaning of each
parameter can guide the backcalculation program during convergence.

Log |E*|

d (in c re a s e )

(a + b )

C (N e g )
C (P o s )
a

Log

R e d u c e d F re q u e n c y

Figure 5: Physical meaning of each parameter on the sigmoidal function
(Kim, 2008)

Relaxation modulus E(t) is another fundamental material property that describes the
viscoelastic behavior of AC. E(t) can be computed directly from |E*| mastercurve using an
interconversion technique as described in the next subsection. It can also be characterized using a
sigmoidal function with 4 parameters, as shown in Eq. (2). The relaxation modulus for the same
mixture above (PG 64-28) at the same reference temperature is shown in Figure 6. It shows that
the parameter group (a, b, c, d) determines the material properties of the AC mixture.

log(Et ) = a +

b
1+

ec + d *log (t R )

where: tR is the reduced time at the reference temperature.
12

1

(2)

Relaxiation Modulus E(t)

Relaxation Modulus (kPa)

1.0E+08
Experiment
Simulated
1.0E+07
1.0E+06
1.0E+05
1.0E+04
1.0E+03
1.E-14 1.E-11 1.E-08 1.E-05 1.E-02 1.E+01 1.E+04 1.E+07 1.E+10 1.E+13 1.E+16
Reduced time (second)

Figure 6: Relaxation modulus |E(t)| of one asphalt concrete mixture (PG64-28)
(Kutay, 2008)

2.3 Interconversion between E(t) and |E*| for AC Mixture
Mastercurve (|E*|), relaxation modulus E(t), and creep compliance D(t) are equivalent for
viscoelastic materials, i.e., there is a unique relationship between |E*|, E(t) and D(t). The
interconversion between |E*| and E(t) is discussed in this section, and the interconversion
between E(t) and D(t) will be discussed in the following section.
In viscoelasticity theory, E(t) and D(t) are calculated from experimental results of
mastercurve |E*|, and are often expressed as the Prony series. It is important to note that a single
polynomial model cannot be used for fitting the entire mastercurve, because the polynomial
swing at low and high temperatures causes irrational modulus value predictions when
13

extrapolating outside the range of data (Kim, 2008). Typically, E(t) curve is represented using the
generalized Maxwell model, shown in Figure 7.

E
o
E
1
E
2

?2

E
n

s

?1

?n

s

Figure 7: The generalized Maxwell model
E*(ω) is mathematically expressed as the storage modulus (E’ (ω)) and loss modulus (E”
(ω)). The relationship between them is given in Eq. (3):
2
N
⎧
( ωTi )
⎪ E ' ( ω) = E * ( ω) *cos ( φ ) = E 0 + ∑ E i
2
1 + ( ωTi )
i =1
⎪
⎨
N
ωTi
⎪E" ω = E * ω *sin φ = E
( )
( ) ∑ i
2
⎪ ( )
1 + ( ωTi )
i =1
⎩

where:

Ti =

(3)

ηi
Ei

, also called relaxation time.

If all the coefficients of Ei and Ti are known, E*(ω) can be directly calculated by Eq. (3).
On the other hand, if E*(ω) and φ(ω) are given, then Ti will be assumed for each component of
14

-15

the general Maxwell model, typically ranging from 10

15

to 10

seconds, and a group of linear

equations is set up to solve for Ei. Usually there are more equations than unknowns, and the
technology of least square root of error (LSRE) may be used to calculate Ei. If the coefficients of
Ei are found, the relaxation modulus is mathematically expressed as:

N

E (t ) = E0 + ∑ Ei e

−

t
Ti

(4)

i =1

2.4 Interconversion between E(t) and D(t) for AC Mixture
The mathematical relationship between E(t) and D(t) in the time domain is given by the
following integral (Kim, 2008):

t

∫ E (t − τ )
0

dD(τ )
dτ ≡ 1
dτ

(5)

Solving the integral for interconversion can be done using a numerical approach. This requires
that the integral be divided into a large number of time segments, which is consistent with the
relaxation modulus from the generalized Maxwell model. This can be done by using the Prony
series form of the E(t) and D(t) functions. Given

{τ

j

{ρi , Ei (i = 1,2,..., m ) and E0 } or

, D j ( j = 1,2,..., m ) and D0 } and the target time constants, the unknown set of

constants can be determined through a system of linear algebraic equations. For example, creep
compliance in its Prony series form, Dj (j=1,2,…,n), can be determined from the relaxation
15

modulus E(t) as follows (Kim, 2008):

[A]{D} = {B}
Or:

(6)

Akj D j = Bk (j=1,2,…,n; k=1,2,…,p) where:

t
⎧ ⎛
− k ⎞
⎪E0 ⎜1 − e τj ⎟ +
⎟
⎪ ⎜
⎝
⎠
⎪
⎪ m ρ E ⎛ − tk − tτk ⎞
⎪∑ i i ⎜ e ρi − e j ⎟ ( ρi ≠ τ j )
⎜
⎟
⎪ i=1 ρi − τ j ⎝
⎪
⎠
Akj = ⎨
t
− k ⎞
⎪ ⎛
τ
E0 ⎜1 − e j ⎟ +
⎪
⎟
⎪ ⎜
⎝
⎠
⎪
t
m
⎛ − ρk ⎞
⎪ t k Ei
⎜ e i ⎟ ( ρi = τ j )
⎪∑
⎜
⎟
i =1 τ j ⎝
⎪
⎠
⎩

16

(7)

m

Bk = 1 −

E0 + ∑ Ei e
i =1

−

tk

ρi

m

E0 + ∑ Ei
i =1

The symbol tk (k=1,2,…,p) represents a discrete time corresponding to the upper limit of
integration in Eq. (5). Once the model constants Dj, and τj are found, the function D(t) can be
obtained in its Prony series form. Similarly, E(t) can be calculated from D(t) by solving for the
unknown constants of the Prony series representation of E(t).
An example of such calculation is given here for an AC mixture. The source and target
Prony series coefficients of a specific AC mixture are shown in Table 1, and the corresponding
fitted data are plotted in Figure 8.

17

Table 1. Prony series coefficients for E(t) and D(t) for a typical AC mixture

Relaxation Modulus
(MPa)
E∞
ρi(s)
1.E-10
1.E-09
1.E-08
1.E-07
1.E-06
1.E-05
1.E-04
1.E-03
1.E-02
1.E-01
1.E+00
1.E+01
1.E+02
1.E+03
1.E+04
1.E+05
1.E+06
1.E+07
1.E+08
1.E+09
1.E+10
1.E+11
1.E+12
1.E+13

Creep Compliance
(1/MPa)

D0
9.8
Prony Coefficients
Ei (MPa)
1.14E+02
6.81E+01
1.81E+02
3.67E+02
6.27E+02
1.02E+03
1.92E+03
2.92E+03
3.82E+03
3.76E+03
2.38E+03
1.28E+03
4.93E+02
1.66E+02
5.59E+01
2.08E+01
8.64E+00
4.11E+00
2.13E+00
9.98E-01
8.00E-01
9.92E-02
3.25E-01
6.99E-02

τi(S)
1.E-10
1.E-09
1.E-08
1.E-07
1.E-06
1.E-05
1.E-04
1.E-03
1.E-02
1.E-01
1.E+00
1.E+01
1.E+02
1.E+03
1.E+04
1.E+05
1.E+06
1.E+07
1.E+08
1.E+09
1.E+10
1.E+11
1.E+12
1.E+13

18

5.20E-08
Di (1/MPa)
3.87E-11
3.03E-10
5.85E-10
7.52E-10
1.97E-09
3.42E-09
6.75E-09
1.35E-08
2.82E-08
6.47E-08
1.48E-07
3.48E-07
1.03E-06
3.02E-06
8.43E-06
1.63E-05
2.05E-05
1.76E-05
1.27E-05
8.80E-06
4.90E-06
3.64E-06
1.27E-06
2.11E-06

Relaxation Modulus (MPa)

E(t) (MPa)
D(t) (MPa)

1.E+04

1.E-04
8.E-05

1.E+03
6.E-05
1.E+02
4.E-05
1.E+01

2.E-05

Creep Compliance (1/MPa)

1.E-04

1.E+05

1.E+00
0.E+00
1.E-15 1.E-11 1.E-07 1.E-03 1.E+01 1.E+05 1.E+09 1.E+13 1.E+17
Reference Time (s)
(S)

Figure 8: E(t) and D(t) from interconversion for a typical AC mixture

2.5 Motivation for the Study
Current backcalculation methods, both static and dynamic methods, ignore the viscoelastic
behavior of AC layer. Assuming AC as a linear elastic material is no longer a valid and
appropriate assumption. The new Mechanistic Empirical Pavement Design Guide (ME-PDG)
utilizes the |E*| mastercurve of AC layer. Therefore, there is a growing need for estimating entire
|E*| mastercurve of AC.
There are three major challenges for the VE backcalculation. First, the dynamic response of
the pavement under field FWD test should be estimated correctly. Second, the elastic modulus of
the base and subgrade layer should be as accurate as possible. Third, the AC material relaxation

19

modulus E(t), or the dynamic modulus mastercurve |E*| should be expressed as a function
instead of one constant value.
The pavement response under the FWD test is the response of a dynamic system after load
impulse, and the phenomenon of wave propagation, as evidenced from the time delay of each
sensor. Additionally, backcalculation algorithms are typically more sensitive to the elastic
modulus of the base and subgrade layer than the AC layer. This indicates that the accuracy of the
modulus of base and subgrade layer significantly affect the accuracy of the AC layer. Based on
the wave propagation theory, the elastic modulus of the semi-infinite subgrade layer (Esubgrade)
can be calculated from the time delay of each sensor, if the interface of different physical layer is
insignificant for wave propagation. This forward calculated modulus can be used as the seed
value of Esubgrade in the VE backcalculation program, and it accelerates the convergence of
backcalculation. A detailed discussion will be given in Chapter 3.
In AC mixture modeling, the generalized Maxwell model is often employed to illustrate the
viscous property of the material. It is the constitutive equation in the time domain. In order to
capture the material constitutive behavior, usually there are more than 30 groups of Maxwell
model as a whole in the physical model. If the backcalculation is done directly to calibrate the
generalized Maxwell model, there would be more than 30 variables in the model, and it is almost
impossible for backcalculation. If the sigmoidal function is chosen, the backcalculation is used to
calculate the parameters of the sigmoidal function (a, b, c, d) in either E(t) or |E*|, although the
VE forward calculation is done in time domain. The detailed explanation of the VE forward
solution is presented in Chapter 4, and the backcalculation algorithm is given in Chapter 5.
The objective of this study is to propose a methodology that is robust in the backcalculation

20

of both AC viscoelastic modulus and the elastic modulus of base and subgrade layer, based on
viscoelasticity theory. More specifically, the thickness, the Poisson’s Ratio and mass density of
each layer are accurately given as input in backcalculation, only the four parameters (a, b, c, d)
of the sigmodial function of E(t) or |E*|, and the elastic modulus of the base (Ebase) and
subgrade (Esubgrade) layers need to be identified.

21

Chapter 3. Estimating the Subgrade Elastic
Modulus using Wave Propagation Method
3.1 Introduction
The Spectral Analysis of Surface Waves (SASW) method is a dynamic non-destructive test
method for determining the shear wave velocity and shear modulus of soils in-situ. The test
procedure was developed in the 1980s at the University of Texas, Austin, and included three
steps: data acquisition, dispersion analysis, and inversion. The test is done at different
frequencies, and the shear modulus is calculated by the frequency and the phase angle difference
between sensors at different locations, through the dispersion curve, which is a plot of phase
velocity vs. frequency. A detailed procedure of this test can be found in Ganji et al. (1998).
A similar method to estimate the elastic modulus of the subgrade layer is proposed in this
chapter. It is based on the spectral element method and the wave propagation theory. The steps of
the proposed procedure are as follows:
Assume reasonable material properties for the subgrade layer, including mass density
and Poisson’s ratio;
Obtain the wave propagation velocity, or Rayleigh wave velocity (VR), based on the
location and the occurrence time of the peak deflection for each sensor in the FWD test;
Obtain the k-value, or the ratio of Rayleigh wave velocity (VR) and shear wave velocity
(Vs);
22

Calculate the shear modulus of the subgrade layer, using Vs and estimated mass
density;
Calculate the elastic modulus, based on shear modulus and Poisson’s ratio.
The mass density is often estimated by the level of compaction during pavement
construction, based on the engineer’s experience; while the Poisson’s ratio is often estimated by
the type of the soil material, i.e., clay, sand, or bedrock. The VR is calculated by the distance
between sensors divided by the difference of the occurrence time of the peak deflection for each
sensor, from the time-history in the FWD test. The only challenge is to find the ratio between VR
and Vs, or the k-value.
For a plane wave expressed in Cartesian coordinates, the ratio (K) of the velocity of
Rayleigh Wave (VR) to the velocity of Shear Wave (VS) is only a function of Poisson’s ratio. K is
the root of the following equation (Richart et al., 1970):

(

)

(

)

K 6 − 8K 4 + 24 − 16α 2 K 2 + 16 α 2 − 1 = 0 ,
where:

α2 =

(8)

1 − 2ν
2 − 2ν

The FWD test generates axial-symmetric waves instead of plane wave, above result may not
be valid, and a proposed method to estimate k-value is briefly discussed in the following section.

3.2 Proposed Method to Estimate k-value
First, the theoretical formulation for the spectral element method is briefly discussed.
Second, the dominant frequency in the FWD test is briefly introduced. Third, numerical
examples are presented to calculate the k-value for specific subgrade site at specific frequencies.
23

Fourth, sensitivity analysis is presented for k-value for a range of elastic moduli, with different
Poisson’s ratios of the subgrade layer, under cyclic loading of different frequencies.

3.2.1 Theoretical Formulation for Spectral Element Method
There are two types of waves in an infinite, homogeneous, isotropic elastic medium: waves
of dilatation and waves of distortion. In a half-space medium, a third wave, the Rayleigh wave, is
found as a solution for the equations of motion. It corresponds to a wave whose motion is
confined to a zone near the boundary of the half-space (Richart et al., 1970).
In a cylindrical coordinate system, assuming a wave traveling in the radial direction, particle
displacement will be independent of the angular direction. The equations of motion of an
isotropic linear elastic material can be expressed in terms of the displacements by use of Navier’s
equations (Al-Khoury, 2001), as:

..

( λ + μ ) ∇∇.u + μ∇ 2 u = ρu ,

(9)

The vector u represents the displacement components of the material, and ρ is the mass
density of the material. ∇ indicates a vector differential operator, ∇ . u is the divergence of u,
and

∇ 2 . u is the Laplace operator of u. For axial symmetry, ∇ 2

∂2 1 ∂ ∂2
+
∇ = 2+
r ∂r ∂z 2
∂r
2

is expressed as:

(10)

where λ and μ are the Lame Constants expressed as:

λ=

νE
,
1 + ν )(1 − 2ν )
(
24

(11)

μ=G=

E
,
2 (1 + ν )

(12)

where E is Young’s modulus and ν is Poisson’s ratio.
In the Helmholtz decomposition, the displacement vector u is expressed as the sum of the
gradients of a scalar potential φ and the curl of a vector potential ψ as:

u= ∇ φ+ ∇ × ψ

(13)

In view of the axial symmetry, the vector ψ only has a component in the θ-direction, which
is

ψ=ψ eθ .

Denoting the displacement components in the r and z directions by u and w,

respectively, the relations between the displacement components and the potentials are shown
below (Yeh, 1989):

∂ϕ ∂ψ
⎧
u=
−
⎪
∂r ∂z
⎪
⎨
⎪ w = ∂ϕ + 1 ∂ ( rψ )
⎪
∂z r ∂r
⎩

(14)

Because of Eq. (9) and (13), the above potentials satisfy the following axial symmetric wave
equations:

⎧ ∂ 2 ϕ 1 ∂ϕ ∂ 2 ϕ 1 ∂ 2 ϕ
+ 2 = 2 2
⎪ 2 +
r ∂r ∂z
cp ∂t
⎪ ∂r
,
⎨ 2
2
2
⎪ ∂ ψ + 1 ∂ψ + ∂ ψ − ψ = 1 ∂ ψ
⎪ ∂r 2 r ∂r ∂z 2 r 2 c2 ∂t 2
s
⎩
25

(15)

where the constants cp and cs are defined as:

⎧
⎪cp =
⎪
⎨
⎪c =
⎪ s
⎩

λ + 2μ
ρ
μ
ρ

(16)

Based on Eq. (15), there are two types of waves: compression wave (P) and Shear wave (S)
respectively. Because of isotropy, these waves are uncoupled from each other.
Considering the harmonic loading case with frequency ω, the potentials can be solved by the
method of separation of variables with some arbitrary constant k, assuming the format of the
potentials are (Maurice et al., 1957):

⎧ ϕ = S ( kr ) T ( z )
⎪
,
⎨
⎪ψ = S1 ( kr ) T1 ( z )
⎩

(17)

provided that the function S and T satisfy the equations:

⎧ d 2S 1 dS ⎛ ω2
⎞
2
+ ⎜ 2 − υ ⎟S = 0
⎪ 2+
⎟
r dr ⎜ c p
⎪ dr
⎝
⎠
,
⎨
⎪ d 2T
2 2
⎪ 2 +υ T =0
⎩ dz

26

(18)

υ =
2

where

ω2
c

2
p

2
− k 2 = kα − k 2

υ' =
with
2

Similar equations hold for S1, T1,

equation of Eq. (18) has a particular solution:

ω2
cs2

2
− k 2 = kβ − k 2 .

The second

T = e − ivz , which remains finite as z -> +∞ .

The first equation of Eq. (18) is a form of Bessel’s equation, which has the Bessel function
solution J0(kr), which remains a non-zero value at r = 0. Thus, two particular solutions are
obtained:

⎧ ϕ = Ae − ivz J 0 ( kr )
⎪
,
⎨
− iv 'z
⎪ψ = Be J1 ( kr )
⎩

(19)

A and B are constants determined by the boundary conditions. The solutions vanish as z—>
∞ and also vanish as r—> ∞, because of the property of the Bessel function J0(kr). The
displacements can be expressed as:

(
(

)
)

⎧u = − Ake− ivz + iv' Be −iv 'z J1 (kr )
⎨
w = − iAve−ivz + Bke−iv 'z J 0 (kr )
⎩
The stress equations in cylindrical coordinates are:

27

(20)

⎧
⎛ ∂u ∂w ⎞
τzr = μ ⎜ +
⎟
⎪
∂z ∂r ⎠
⎪
⎝
⎨
⎪σ = λ + 2μ ∂w + λ ∂ ( ru )
)
⎪ zz (
∂z r ∂r
⎩

(21)

By replacing u and w, the following equation is obtained at the surface z = 0:

(

)

2
⎧ τrz = ⎡ 2ikυAu + k β − 2k 2 uB⎤ J1 ( kr )
⎪
⎣
⎦
⎨
2
2
⎡
⎤
⎪σzz = ⎣ − k β − 2k uA − 2ikυ 'uB⎦ J 0 ( kr )
⎩

(

)

(22)

To make the problem discrete, some boundary conditions in the radial direction are
introduced. At the radial boundary r = R (far away from the source), the amplification of the
oscillation is considered to vanish. Only the vertical displacement w is important in the FWD test,
so the horizontal component will be ignored in the future. These considerations can be
implemented by the infinitely many positive roots

αm

of the

J0

function, as

km= αm/R

(Al-Koury et al., 2001). For FWD test loading, if the maximum pressure is P, the appropriate
boundary conditions are (Maurice et al., 1957):

⎧ τrz z =0 = 0
⎪
⎨
⎪σzz z =0 = PJ 0 ( kr )
⎩
By combining Eq. (22) and Eq. (23), the coefficients A and B can be solved.

28

(23)

(

)

2
⎧
k β − 2k 2 P
⎪A =
F(k) μ
⎪
,
⎨
⎪ B = − 2ikυ ' P
⎪
F(k) μ
⎩

(

F (k ) = 2k − k β
where:
2

)

2 2

(24)

− 4k 2υυ ' , which is called Rayleigh’s function

(Maurice et al, 1957). The expression of the wave propagation includes Rayleigh’s function, so
this wave is called Rayleigh wave (Xiang, 2009).
Fig. 9 shows different kinds of wave propagation on the surface of an ideal medium after a
load impulse. It illustrates the vibration magnitude by different wave modes, and there is a tiny
vibration for compression wave (P) just after the load impulse. However, the vibration of shear
wave (S) and Rayleigh wave (R) happen almost simultaneously, especially if the location is near
to the loading point. The time difference of S wave and R wave is insignificant as the velocities
of S wave and R wave are almost identical. Thus, the distance between two adjust peaks in
vertical displacement is considered the wave length. However, wave length varies by location or
propagation time, and the velocity is hard to identify, because the interaction of shear wave and
Rayleigh wave affects the peak displacement.

29

t

t

Figure 9: Wave system from surface point source in ideal medium
(a) Effect of different wave modes on horizontal displacement
(b) Effect of different wave modes on vertical displacement
1
(c) The combined vibration for particle at Point ○
(Richart et al., 1970)

Typically, the S wave or R wave velocity is around 150 m/s (500 ft/s) for subgrade layers,
and the effective frequency in FWD test is less than 50 Hz, i.e., the S wave length is more than 3
m (10 ft), which is beyond the range of the sensor distribution in the FWD test. Thus, the
occurrence time of peak deflection for each sensor is considered, noted as the time delay, rather
than the deflection basin, and by the difference of the occurrence time, the velocity of R wave is
identified.
30

3.2.2 Dominant Frequency of the FWD Test
Eq. (24) indicates the Rayleigh wave propagation velocity is a function of frequency, while
the shear wave velocity is independent to frequency. K-value should be dependent on the loading
frequency. Thus, the equivalent frequency of the FWD test, or the dominate frequency range in
the FWD test should be identified.
The equivalent frequency of a loading history can be identified by: (a) the centroid of the
area formed by Fourier Series; (b) the centriod of the area formed by the power spectral density
(PSD) in the frequency domain; (c) the central frequency method from earthquake engineering,
as defined below. In method (a), the coefficient of Fourier series is found first, and then the
centroid of the coefficient corresponding to the frequency is obtained as the equivalent frequency.
Method (b) is well-known as Parseval’s theorem or the principle of energy conservation in time
and frequency domains. The power spectral density function (G (ω)) is calculated by coefficients
of Fourier series first, and the centroid of the PSD is called the equivalent frequency. In method

(c), the equivalent frequency is calculated as:

ω=

λ2
λ0

where λi

=∫

ωN

0

w i G (ω )dω .

For a typical FWD loading history, the loading duration is 35 ms with a haversine function.
It is assumed that the rest time is 1 s with a sampling interval of 1ms, which only affects
accuracy by a numerical error of 2 - 3 %. The equivalent frequency is calculated by these
methods as: (a) f = 17.23 Hz. (b) f = 13.25Hz. (c) f = 16.38 Hz. Considering the actual loading is
not a haversine function, and the loading time may differ from 35 ms, the equivalent frequencies
for the three methods are listed in Table. 2. It shows the typical range of the dominant frequency
in the FWD test is 5 - 25 Hz.

31

Table 2: The equivalent frequency for different load duration times

Duration time (ms)
25
30
35
40
45
50
Method (a)
24.09 Hz 20.10 Hz 17.23 Hz 15.09 Hz 13.42 Hz 12.08 Hz
Method (b)
18.56 Hz 15.47 Hz 13.25 Hz 11.60 Hz 10.31 Hz 9.28 Hz
Method (c)
22.93 Hz 19.11 Hz 16.38 Hz 14.33 Hz 12.74 Hz 11.46 Hz

3.2.3 Numerical Example to Calculate k-value
Eq. (24) indicates the VR under cyclic loading depends on the loading frequency, elastic
modulus and Poisson’s ratio. One example of wave propagation is illustrated for the hypothetical
profile which only includes the half-space subgrade layer. The cyclic loading frequency is 5 HZ,
and the material properties are: elastic modulus = 68.95 MPa, and Poisson’s ratio = 0.45. The
phase angle of each point at the surface is shown in Fig. 10 for t = 0 second.
Phase angle at the surface of the subgrade layer at t=0
(E = 68.95 Mpa, v = 0.45, f = 5 Hz)

Phase angle (o)

90
60
30
0
-30
-60
-90
0

10

20
30
Location (m)

40

50

Figure 10: Phase angle of each point at the surface for a specific example (t =0 s)

The wavelength is defined as the spatial period of the phase angle, or the distance over
32

which the phase angle’s shape repeats. The wavelength in the above example is around 18.5 m,
and the wave velocity (VR) is around 92.5 m/s, as the loading frequency is 5 Hz.
The mass density (ρ) of the subgrade layer is assumed as 2000 kg/m3, and the shear modulus
(G) is 23.78 MPa, calculated from elastic modulus and Poisson’s ratio (Eq. 12). The shear wave
velocity is 109 m/s, from Eq. (25).

G
ρ

Vs =

(25)

The k-value is 0.8486, which is VR divided by Vs (Eq. 26).

k=

VR
Vs

(26)

The k-value only depends on the Poisson’s ratio in plane wave, independent of loading
frequency and elastic modulus of the subgrade layer. Similar conclusion is expected for the axial
symmetric FWD test.
Cyclic loading is applied in the SASW method; however, only load impulse is applied in the
FWD test. The phase angle for each point at the surface of the subgrade layer can be replaced by
the occurrence time of the peak deflection.

3.2.4 Sensitivity Analysis for k-value
The typical range of the elastic modulus of subgrade is 35 - 280 MPa (5 to 40 ksi) (NCHRP,
2004). The dominant frequency in the FWD test is 5 to 25 Hz. Since the typical Poisson’s ratio
ranges from 0.20 to 0.49, a sensitivity analysis for k-value with different combinations of elastic

33

modulus and frequency was conducted for possible values of Poisson’s ratio (Fig. 11 - 17).

R

S

K−value (V /V )

1
0.8
0.6
0.4
0.2
0
30
8

20

6
4

10
0

Frequency (Hz)

8

2

x 10

0
Elastic Modulus (MPa)

Figure 11: Sensitivity analysis for k-value (ν= 0.20)

0.8

R

S

k−value (V /V )

1

0.6
0.4
0.2
0
30
8

20

6
4

10
Frequency (Hz)

2
0

0
Elastic Modulus (MPa)

Figure 12: Sensitivity analysis for k-value (ν= 0.25)

34

8

x 10

0.8

R

S

k−value (V /V )

1

0.6
0.4
0.2
0
30
8

20

6
4

10
0

Frequency (Hz)

8

2

x 10

0
Elastic Modulus (MPa)

Figure 13: Sensitivity analysis for k-value (ν= 0.30)

0.8

R

S

k−value (V /V )

1

0.6
0.4
0.2
0
30
8

20

6
4

10
Frequency (Hz)

2
0

0
Elastic Modulus (MPa)

Figure 14: Sensitivity analysis for k-value (ν= 0.35)

35

8

x 10

k−value (VR/VS)

1
0.8
0.6
0.4
0.2
0
30
8

20

6
4

10
0

Frequency (Hz)

8

2

x 10

0
Elastic Modulus (MPa)

Figure 15: Sensitivity analysis for k-value (ν= 0.40)

k−value (VR/VS)

1
0.8
0.6
0.4
0.2
0
30
8

20

6
4

10
Frequency (Hz)

2
0

0
Elastic Modulus (MPa)

Figure 16: Sensitivity analysis for k-value (ν= 0.45)

36

8

x 10

k−value (VR/VS)

1
0.8
0.6
0.4
0.2
0
30
8

20

6
4

10
0

Frequency (Hz)

8

2

x 10

0
Elastic Modulus (MPa)

Figure 17: Sensitivity analysis for k-value (ν= 0.49)

Sensitivity analysis indicates the k-value is almost constant if the elastic modulus of the
subgrade is less than 500 MPa (70 ksi); it is not sensitive to the frequency of the loading for the
dominant range in the FWD test, but it is slightly affected by Poisson’s ratio. Thus, one average
k-value could be used for each Poisson’s ratio. Table 3 shows the average k-value of all the entire
sensitivity range as well as for the more practical range of moduli (5 to 40 ksi) and frequencies (5
to 25 Hz).
Table 3: The average k-value for different cases of Poisson’s ratio

Poisson’s ratio
Average
Entire range
k-value Practical range
COV
Entire range
(%)
Practical range

0.20
0.7957
0.7955
1.63
0.83

0.25
0.8051
0.8056
1.58
0.76

0.30
0.8132
0.8140
1.65
0.94

37

0.35
0.8176
0.8189
1.89
1.22

0.40
0.8173
0.8197
2.34
1.78

0.45
0.8102
0.8154
3.32
2.77

0.49
0.7970
0.8073
5.48
4.66

Table 3 shows the average k-value does not significantly affected by Poisson’s ratio. It
would be much better to use a single average k-value to simplify the problem; k=0.8109 is
suggested. However, it should be mentioned that the k-value may be quite different if there is a
stiff layer underneath, as the cases are not verified for subgrade with elastic modulus higher than
700 MPa (100 ksi).

3.3 Proposed Procedure to Estimate Subgrade Elastic Modulus
Based on the result in Section 3.2, the following steps are applied to estimate the elastic
modulus of the subgrade layer using k-value.
1. Estimate the Poisson’s ratio (ν) and mass density (ρ) of the subgrade layer. The
Poisson’s ratio depends on the soil type, and it is well documented in the literature. The
soil type, compaction level during construction, and water table level will affect the
3

mass density, ranging from 1500 to 2000 kg/m . However, the estimated elastic
modulus is not significantly affected by the mass density, and the engineers’ experience
will be helpful for the estimation.
2. Find the slope in the plot of time delay vs. location. Six sensors are located at different
positions in the FWD test, and the occurrence time of peak deflection for each sensor is
different. Time delay of each sensor is defined as the difference between the time of
peak loading and the time of peak deflection. Time delay and locations can be plotted in
one figure, and the slope (m) for the farther sensors is found by linear regression.
One example of the load impulse and time-history for the FWD test are shown in Fig.
18. The test was conducted at LTPP (Long Term Pavement Performance) station 04-1036 in
1998. The plot of time delay vs. location is shown in Fig. 19. The slope in Fig.19 is
38

calculated as: m = 0.0069 s/m.

FWD test history

Load (kPa)/deflection (um)

800
load
Sensor
Sensor
Sensor
Sensor
Sensor
Sensor
Sensor

700
600
500
400
300

1
2
3
4
5
6
7

200
100
0

-100

0

10

20

30
Time (ms)

40

50

60

Figure 18: Example of FWD load and deflection time histories at LTPP station 04-1036 in
1998

39

Time delay vs. Location
0.01
0.009
Time delay (ms)

0.008
0.007
0.006
y = 0.0069x - 0.0013
2
R = 0.9978

0.005
0.004
0.003
0.002
0.001
0
0.0

0.2

0.4

0.6

0.8
1.0
Location (m)

1.2

1.4

1.6

1.8

Figure 19: Time delay vs. location for the example of FWD test at LTPP station 04-1036 in
1998

3. Calculate the Rayleigh wave velocity (VR). The wave propagation velocity is the
inverse of the slope (m) in the plot of time delay vs. location. In the current example, VR
= 1/0.0069 = 144.9 m/s.
4. Calculate the Shear wave velocity (Vs). The Vs is VR divided by k-value. In this
example, Vs = 144.9/0.8109 = 178.1 m/s.
5. Calculate the shear modulus (G) of the subgrade layer. The mass density is estimated in
step 1, as ρ=1800 kg/m3. The relationship between G and Vs is shown in Eq. (27),
which is similar to Eq. (25) in Section 3.2.

40

G = Vs2 * ρ

(27)

The shear modulus is calculated as: G= 57.1 MPa in the current example.
6. Calculate the elastic modulus (E) of the subgrade layer by Eq. (28), with assumed
Poisson’s ratio (ν) in Step 1.

E = 2 (1 + υ ) G

(28)

The elastic modulus in this example is E = 165.6 MPa, assuming ν = 0.45.
For simplification, a flow chart for this procedure is presented in Fig. 20.
Assume ρ and ν
Find time delay and location for each sensor in FWD
Plot time delay vs. location, as Fig. 16
Obtain the slope (m) for farther sensors

Calculate Rayleigh wave velocity, VR=1/m
Calculate shear wave velocity, VS=VR/k (k=0.8109)
Calculate shear modulus, by Eq. (27)
Calculate elastic modulus, by Eq. (28)
Figure 20: Flow chart to estimate subgrade elastic modulus using k-value

It is difficult to identify how many sensors should be used in step 2. Often there are 6 - 7
41

sensors in the FWD test, so the farther sensors would be sensors 4 - 6 or 5 - 7. There are
interactions between interfaces of different layers, so the time delay for close sensors is
influenced by the combined effect of different layers. Farther sensors are less affected by
interactions, and it is preferable to obtain the m-value from multiple sensors. The slope between
the farthest sensors is the dominant factor for m-value, and the time delay for other sensors
2

should be ignored if including them reduces the regression R to below 0.98.
The proposed procedure uses the time delay of each sensor in the FWD test, while the
SASW method uses the phase angle of each sensor. Although load impulse can be converted to
summation of different magnitudes of cyclic loading at different frequencies by Fourier
Transform, the phase angle for each frequency is very sensitive in converting. Furthermore,
FWD time histories are typically truncated before they fully decay, leaving errors when using
field measurements. The pavement response in the frequency domain is not as accurate as in the
time domain. Therefore, the time delay instead of phase angle is used in the proposed method.

3.4 Verification of the Proposed Method
No field data are available with all the information on the pavement profile, so only
numerical examples are used to validate the proposed method. The time-history of the field FWD
test is simulated by two well-known dynamic solutions: SAPSI and LAMDA. Both are briefly
introduced first, and numerical examples of subgrade layer with and without ground water table
are presented later, to illustrate the theoretical error of the proposed method.

3.4.1 SAPSI
The forward program SAPSI (Chen, 1987) models the pavement structure as a system of
42

finite layers that are infinite in the horizontal direction and underlain by an elastic half-space
with viscous boundaries. The finite layer solution is based on Kausel’s formulation (Kausel, 1981;
Doyle, 1997) which subdivides the medium into discrete layers that have a linear displacement
function in the vertical direction (finite element method with lumped mass formulation) and
satisfy the wave equation in the horizontal direction (exact formulation). The solution is based on
the premise that if the sublayer thickness is small relative to the wavelength of interest, it is
possible to linearize the transcendental functions and reduce them to algebraic expressions. The
materials are assumed to be isotropic and linearly elastic with hysteretic damping. Full interface
bonding is assumed at the layer interfaces. The mass densities and elastic moduli change with
depth, from layer to layer, but are assumed to be constant within each layer. The top layer
represents the asphalt concrete surface, which can be modeled as a viscoelastic material by
allowing its (complex) modulus to be a function of frequency.

3.4.2 LAMDA
Al-Khoury et al. (2001) developed an efficient axial-symmetric forward solution, called
LAMDA, for the dynamic analysis of flexible pavements using the spectral element method for
the simulation of wave propagation in layered systems. The spectral element method developed
by Doyle (1997) combines elegantly the exact solution of wave motions with the finite element
organization of the system matrices. The system is solved by double summation over the
involved frequencies and the wave numbers (Rizzi and Doyle, 1992). The double summation
approach using Fourier series is computationally advantageous over Kausel’s formulation, which
relies on the numerical evaluation of integrals between zero and infinity. This type of integration
involves singularities if the system has no damping or very sharp peaks for small damping, and it
requires considerable computational time and capacity. The mass distribution is modeled exactly
43

and hence only one element is sufficient to describe a whole layer without the need for
subdivisions. This makes the resulting system of dynamic equations very small and hence
computationally efficient.

3.4.3 Case of Subgrade without Ground Water Table
One specific numerical example of a pavement without ground water table (GWT) subjected
to the FWD test is analyzed in SAPSI, and simulated by the Finite Element Software ABAQUS.
Table 4 presents a summary of the pavement parameters.
Table 4: Basic information of a three-layer pavement structure without GWT

Physical Layer
AC
Base
Subgrade

Elastic Modulus
(MPa (ksi))
Experimental data
150 (21.8)
100 (14.5)

Poisson’s
ratio
0.35
0.35
0.45

Mass density
3
(kg/m )
2300
2000
1500

Thickness
(m (in))
0.1 (4)
0.3 (12)
Infinity

The FWD load is assumed as a haversine function with duration time of 50 ms, and the
maximum pressure is 707 kPa with a loading plate radius of 0.15 m. In SAPSI, 3% damping ratio
is applied for the subgrade layer, as it is a typical value for clay and can eliminate the free
vibration of the system after dynamic loading. The time delay of each sensor is given in Table 5.
Table 5: Time delay of each sensor simulated by SAPSI

Sensor
1
Physical distance (m) 0.00
Time delay (s)
0.0039

2
0.41
0.0049

3
0.67
0.0067

4
0.93
0.0089

5
1.19
0.011

6
1.44
0.0129

7
1.70
0.0148

Since the sampling interval is 0.1 ms in the numerical simulation, the error from calculation
may be relatively high. For the elastic modulus of the subgrade layer, the time delay and the
location of sensors 4 - 7 are used given the thickness of the pavement structure is 0.4 m (16 in).

44

The linear regression for the time delay of these sensors is shown in Fig. 21.

Time delay VS. Location
0.016
y = 0.007656 x + 0.001832
2
R = 0.999474

Time delay (s)

0.014
0.012
0.01
0.008
0.006
0.004
0.002
0
0

0.2

0.4

0.6
0.8
1
1.2
Sensor Location (m)

1.4

1.6

1.8

Figure 21: Time delay vs. location for a numerical example of FWD test (without GWT)

From the statistical analysis of these sensors, the slope is obtained directly by SAPSI, as m =
0.007656 s/m. The Rayleigh wave velocity is the inverse of the slope m, i.e., the Rayleigh wave
velocity VR = 1/m = 130.6 m/s.
The Poisson’s ratio of the subgrade layer is assumed as υ = 0.45 , which is a typical value
for clay. The average k-value equals 0.8109. Thus, Vs= VR/k = 160.5 m/s. The mass density of

the subgrade layer is assumed as

ρ = 1500 kg / m3

, so the shear modulus of the

subgrade layer is calculated in Eq. (18) as G= 38.6 MPa. From elasticity theory, the elastic
modulus E is calculated in Eq. (14) as E= 112.1 MPa. Given the input parameters of the
45

pavement structure (Table 3), this value of E is relatively accurate with an error of 12.1%. This
level of error is acceptable in pavement engineering.

3.4.4 Case of Subgrade with GWT
One specific numerical example of a pavement with GWT during the FWD test is analyzed
in SAPSI, with GWT 2 ft below the surface of the subgrade layer. The material properties of the
saturated subgrade layer are calculated by the assumption that the shear wave velocity is the
same as that in the dry subgrade layer, while the compression wave velocity equals the sound
propagation velocity in water (1450 m/s). Table 6 gives a summary of the parameters of the
pavement profile with GWT.
Table 6: Basic information for a three-layer pavement structure with GWT

Elastic Modulus
(MPa (ksi))
AC
Experimental data
Base
108 (15.0)
Dry subgrade
34.5 (5.0)
Saturated subgrade
42.9 (6.2)
Physical Layer

Poisson’s
ratio
0.35
0.30
0.45
0.4979

Mass density
3
(kg/m )
2300
2000
1500
1682

Thickness
(m (in))
0.15 (6)
0.3 (12)
0.6 (24)
Infinity

The FWD load is assumed as a haversine function with duration time of 35 ms, and the
maximum pressure is 707 kPa with a loading plate radius of 0.15 m. In LAMDA, 3% damping
ratio is applied for the subgrade layer, as it is a typical value for clay and can eliminate the free
vibration of the system after dynamic loading. The simulated time-history of the field FWD test
is shown in Fig. 22.

46

Deflection history with GWT (2ft)
6.E-04

sensor 1
sensor 3

Deflection (m)

sensor 6

sensor 7

4.E-04

sensor 4

sensor 5

5.E-04

sensor 2

sensor 8

3.E-04
2.E-04
1.E-04
0.E+00
-1.E-04
0.01

0.02

0.03

0.04

0.05
0.06
time (sec)

0.07

0.08

0.09

Figure 22: LAMDA simulated time-history due to the FWD test (with GWT)

The time delay of each sensor is given in Table 7.
Table 7: Time delay of each sensor by LAMDA

Sensor
1
Distance (m) 0.00
Delay (s)
0.0051

2
0.20
0.0055

3
0.30
0.0058

4
0.46
0.0065

5
0.61
0.0074

6
0.91
0.01

7
1.22
0.0137

8
1.52
0.0175

Since the sampling interval is 0.1 ms in numerical simulation, the error from calculation may
be relatively high. For the elastic modulus of the subgrade layer, the time delay and the location
of sensors 6 - 8 are considered, given the thickness of the pavement structure 0.45 m (18 in). The
linear regression for the time delay for these sensors is plotted in Fig. 23.

47

Time delay vs. location (GWT=2ft)

Time Delay (s)

0.02
y = 1.2303E-02x - 1.2667E-03
0.015

R2 = 9.9994E-01

0.01
0.005
0
0

0.5

1

1.5

2

Location (m)
Figure 23: LAMDA simulated time-delay for the FWD test (with GWT)

From the statistical analysis of these sensors, the slope is obtained directly by LAMDA, as m
= 0.00123 s/m. The Rayleigh wave velocity is the inverse of the slope m, i.e., the Rayleigh wave
velocity VR = 1/m = 81.3 m/s.
The Poisson’s ratio of the subgrade layer is assumed as υ = 0.45 , which is a typical value
for clay. Using the average k-value of 0.8109, we obtain Vs = VR/k = 99.9 m/s. The mass density
of the subgrade layer is assumed as

ρ = 1500 kg / m3

, so the shear modulus of the subgrade

layer is calculated in Eq. (25) as G = 15.0 MPa. From elasticity theory, the elastic modulus E is
calculated in Eq. (11) as E = 43.4 MPa. Given the input parameters of the pavement structure
(Table 6), this value of E is almost exactly the elastic modulus of the saturated subgrade layer,
with an error of 1.0%. Although the estimated modulus is higher than the modulus of the dry
48

subgrade layer, with an error of about 27.6%, this estimated value is commonly accepted as
reasonable in pavement engineering.
These two numerical examples show that the proposed method to estimate elastic modulus
of the subgrade layer is valid in both cases (without and with GWT), and the validation is
checked by field FWD time-histories from the Long-Term Pavement Performance (LTPP)
database, although the actual elastic modulus of the subgrade layer is unknown.

3.5 Validating the Proposed Method Using Case Studies from LTPP
Three cases are chosen from the LTPP database in this section, and their pavement profiles
and time-histories are shown in each sub-section. Two different methods are used to estimate the
modulus of the subgrade layer. One is calculated by the most common backcalculation software,
MODCOMP5, which is based on minimizing the root mean squared error (RMSE) of the
deflection basin. The other is the proposed method using Rayleigh wave velocity, which is based
on the time delay of each sensor. Since the proposed method works very well in numerical
examples, it is expected that the estimation of the elastic modulus from Rayleigh wave is better
than that from the backcalculation software MODCOMP5.

3.5.1 Example of FWD Test at LTPP 04-1036 Station in 1998
One example of a FWD test at 04-1036 station (Arizona) in 1998 was selected for analysis,
and the FWD test load and deflection histories are shown in Fig. 18 in section 3.3. It illustrates
that the deflection histories yields information on viscoelasticity. There are several overlays
above the original pavement; however, they are combined into one layer as they are all thin AC
mixtures. The maximum pressure during the FWD test is 749 kPa (108.63 psi).

49

The time delay and maximum deflection of each sensor for this example are summarized in
Table 8. Note that the measured time delay values are not as accurate because of the sampling
time interval being too coarse (0.4 ms).
Table 8: Time delay and peak deflection of each sensor due to the FWD test at station
04-1036 in 1998

Sensor
Distance (in)
Time delay (ms)
Peak deflection (mils)

1
0.00
1.6
21.69

2
8.0
1.6
17.05

3
12.0
2.4
13.86

4
18.0
2.8
9.13

5
24.0
2.8
6.50

6
36.0
5.2
3.35

7
60.0
9.2
1.65

Using the wave velocity method proposed in this chapter, the Rayleigh wave velocity is
144.2 m/s, and the elastic modulus is 156.9 MPa (22.7 ksi), assuming the mass density of the
3

subgrade layer is 1800 kg/m .
By MODCOMP5, the elastic modulus of each layer can be calculated, given the pavement
structural profile. The pavement structure information and the back-calculated results from
MODCOMP 5 are shown in Table 9.
Table 9: Structural information and MODCOMP5 backcalculated modulus at LTPP
station 04-1036 in 1998

Physical
Thickness
Layer
(m (in))
AC
0.122 (4.8)
Base (GB)
0.472 (18.6)
Subgrade (SS)
Infinity

Poisson’s
ratio
0.35
0.30
0.45

Mass density
3
(kg/m )
2300
2000
1800

Elastic Modulus
(MPa (ksi))
4771 (692)
86.9 (12.6)
222.7 (32.3)

The deflection basins from the FWD test and the backcalculation result by MODCOMP5 are
shown in Fig. 24. It is clear that the deflection basin is simulated reasonably well, and the
deflections of close sensors match almost exactly. Hence, the backcalculation program converges,
although the result may be unrealistic.

50

Deflection Basin

Maximum Deflection (mils)

0
5
ModComp
FWD test

10
15
20
25
0

10

20

30
40
Sensor Location (in)
Radius distance(in)

50

60

70

Figure 24: Deflection basin from FWD test and MODCOMP5 Simulation at LTPP station
04-1036 in 1998

From the LTPP database, the resilient modulus of the subgrade layer is tested at different
confining pressures and different axial loads, and the summary of the test results is given in
Table TST_UG07_SS07_WKSHT_SUM in the Test Module of the LTPP database. Typically, the
confinement of the subgrade layer is 41.3 kPa (6 psi) in design (Federal Highway Administration,
1996), and the corresponding resilient modulus can be calculated as 102 MPa (14.8 ksi) by
interpolation.
The elastic modulus from the backcalculation by MODCOMP5 is 222.7 MPa (32.3 ksi), and
the modulus from the wave velocity method proposed here is 156.9 MPa (22.7 ksi). Thus, the
difference between the experimental data and the backcalculation result is roughly 118.3%, while
the difference between the experimental data and the result from the proposed wave velocity
51

method is 53.4%. These values are summarized in Table 10.
Table 10: Elastic modulus of the subgrade layer at LTPP station 04-1036 in 1998

Lab_ measured
MR @ 6psi
confinement
102 MPa
(14.8 ksi)

Predicted by
k-value
(Ek-value)
156.9 MPa
(22.7 ksi)

Difference
MODCOMP5
Difference
backcalculation between MR and between MR
Ek-value
(EMODCOMP)
and EMODCOMP
222.7 MPa
53.4 %
118.3 %
(32.3 ksi)

3.5.2 Example of FWD Test at LTPP 32-0101 Station in 1996
Another example is chosen from LTPP station 32-0101 in 1996. The FWD test load and
deflection histories are shown in Fig. 25. There are 5 relatively thick physical layers. For the
backcalculation analysis, they are separated into each physical layer first, and then they are
combined together if the results from backcalculation are unrealistic. The test was done in
Nevada in early winter. The maximum test pressure was 454 kPa (65.85 psi).

52

FWD test history

Load (kPa)/deflection (um)

500

load
Sensor
Sensor
Sensor
Sensor
Sensor
Sensor
Sensor

400
300
200

1
2
3
4
5
6
7

100
0
-100 0

10

20

30
Time (ms)

40

50

60

Figure 25: Example of FWD test results at LTPP station 32-0101 in 1996

The time delay and maximum deflection of each sensor at station 32-0101 in 1996 are
summarized in Table 10. Again, it is worth noting that the measured time delays are not as
accurate, because of the large sampling time intervals (0.4 ms).
Table 10: Time delay and maximum deflection of sensors during the FWD test at LTPP
station 32-0101 in 1996

Sensor
Distance (in)
Time delay (ms)
Peak Deflection (mils)

1
0.00
1.2
4.84

2
8.0
1.6
4.06

3
12.0
1.6
3.39

4
18.0
2.0
2.56

5
24.0
2.0
1.89

6
36.0
5.2
1.18

7
60.0
8.8
0.79

Using the proposed method, the Rayleigh wave velocity is 138.5 m/s, and the elastic
modulus is 120.6 MPa (17.5 ksi), assuming the mass density of the subgrade layer is 1500
3

kg/m .
53

The pavement structural information and the back-calculated results from MODCOMP 5 are
shown in Table 11.
Table 11: Basic information and MODCOMP5 backcalculated modulus at LTPP station
32-0101 in 1996 (Trial: I)

Physical
Layer
AC
Base (GB)
Subbase (GS)
Subbase (TS)
Subgrade (SS)

Thickness
(m (in))
0.183 (7.2)
0.216 (8.5)
0.579 (22.8)
0.305 (12)
Infinity

Poisson’s
ratio
0.35
0.30
0.45
0.45
0.45

Mass density
3
(kg/m )
2300
2000
1800
1800
1500

Elastic Modulus
(MPa (ksi))
7997 (1160)
60.7 (8.8)
5405 (784)
250969 (36400)
165 (24.0)

The result of backcalculation for the subbase (TS) seems to be unreasonable. Another trial
was conducted, in which all the base and subbase layers were combined into one single base
layer with a total thickness of 1.10 m (44.3 in). The results are shown in Table 12.
Table 12: Basic information and MODCOMP5 backcalculated modulus at LTPP station
32-0101 in 1996 (Trial: II)

Physical Layer
AC
Base (combined)
Subgrade (SS)

Thickness
(m (in))
0.183 (7.2)
1.10 (43.3)
Infinity

Poisson’s
ratio
0.35
0.40
0.45

Mass density
3
(kg/m )
2300
2000
1500

Elastic Modulus
(MPa (ksi))
6640 (963)
298 (43.2)
281 (40.8)

The modulus is more reasonable in this trial, although the deflection basin may generate a
higher RMSE. The deflection basins from the FWD test and the backcalculation results obtained
by MODCOMP5 are shown in Fig. 26. It illustrates that the deflection basin is simulated
reasonably well.

54

Deflection Basin

Maximum Deflection (mils)

0
1
ModComp-1
FWD test
ModComp-2

2
3
4
5
6
0

10

20

30
40
Sensor distance (in)
Radius Location (in)

50

60

70

Figure 26: Deflection basin due to the FWD test and two cases of MODCOMP5
simulation at LTPP station 32-0101 in 1996

Similar to the previous example, the subgrade resilient modulus is estimated from the
triaxial laboratory test results to be 64.0 MPa (9.3 ksi) at 6 psi confinement. The elastic modulus
obtained from backcalculation by MODCOMP5 is 281 MPa (40.8 ksi), and the elastic modulus
by the wave velocity method is 120.6 MPa (17.5 ksi). Thus, the difference between the
laboratory experimental data and the backcalculation is 339.1%, as compared to 88.2% when
using the wave velocity approach. These values are summarized in Table 13.
Table 13: Elastic modulus of the subgrade layer at LTPP station 32-0101 in 1996

Lab_ measured
MR @ 6psi
confinement
64.0 MPa
(9.3 ksi)

Predicted by
k-value
(Ek-value)
120.6 MPa
(17.5 ksi)

Difference
MODCOMP5
Difference
backcalculation between MR and between MR and
(EMODCOMP)
Ek-value
EMODCOMP
281 MPa
88.2%
339.1%
(40.8 ksi)
55

The difference between the experimental data and the result from proposed wave velocity
method is still high, because the total thickness of the pavement structure is 50.5 inches, and the
distance for the farthest sensor in the FWD test is only 60 inches, so the time delay for sensor 7
does not directly indicate the elastic modulus of the subgrade layer, rather it is the combination
of the subbase and the subgrade layers.
The freezing Index for this station is 246.63 degree-day, and the frost penetration depth is
about 25 inches by interpolation from the chart (Yoder and Witczak, 1975). The thickness of the
total pavement structure is more than 50 inches, so the frost problem should not affect the
subbase (TS) and subgrade layer.
Since only one FWD test is input in MODCOMP5, the nonlinear property of the subgrade
layer cannot be obtained, and the subgrade layer is assumed as linear elastic material, or the
modulus is independent of the stress condition. This simplification would significantly increase
the error from MODCOMP5.

3.5.3 Example of the FWD Test at LTPP Station 06-0565 in 1999
Another example is chosen from LTPP station 06-0565 in California in 1999. The FWD test
load and deflection histories are shown in Fig. 27. There are 5 relatively thin AC layers, so they
are combined into one AC layer for analysis, although they were constructed at different times.
The maximum applied pressure is 532 kPa (77.16 psi).

56

FWD test history

Load (kPa)/deflection (um)

1000

load
Sensor
Sensor
Sensor
Sensor
Sensor
Sensor
Sensor

800
600
400
200

1
2
3
4
5
6
7

0
-200 0

10

20

30
Time (ms)

40

50

60

Figure 27: Example of the FWD test history at LTPP station 06-0565 in 1999

The time delay and maximum deflection of each sensor are summarized in Table 14.
Table 14: Time delay and maximum deflection of the sensors due to the FWD test at
LTPP station 06-0565 in 1999

Sensor
Distance (in)
Time delay (ms)
Peak deflection (mils)

1
0.00
3.6
32.56

2
8.0
4.4
26.18

3
12.0
5.2
20.75

4
18.0
6.4
15.04

5
24.0
8.0
11.02

6
36.0
11.2
7.24

7
60.0
16.0
4.41

Following the proposed method, the Rayleigh wave velocity is 116.0 m/s, and the elastic
modulus is 101.5 MPa (14.7 ksi), assuming the mass density of the subgrade layer is 1800 kg/m3.
The pavement structure information and the back-calculated results from MODCOMP 5
(trial I) are shown in Table 15.

57

Table 15: Structural information and MODCOMP5 backcalculated modulus at LTPP
station 06-0565 in 1999 (Trial: I)

Physical
Thickness
Layer
(m (in))
AC
0.208 (8.2)
Base (TB)
0.119 (4.7)
Subbase (GS) 0.544 (21.4)
Subgrade (SS) Infinity

Poisson’s
ratio
0.35
0.30
0.45
0.45

Mass density
3
(kg/m )
2300
2000
1800
1800

Elastic Modulus
(MPa (ksi))
834.3 (121.0)
6.9 (1.0)
344738 (50000)
9.7 (1.4)

The result of backcalculation for the subbase (GS) is unreasonable, partially due to the soft
subgrade, and partially due to the soft base (TB) layer. Another trial is conduced by combining
all the base and subbase layers into one single base layer. The results for trial II are shown in
Table 16.
Table 16: Structural information and MODCOMP5 backcalculated modulus at LTPP
station 06-0565 in 1999 (Trial: II)

Physical
Layer
AC
Base (combined)
Subgrade (SS)

Thickness
(m (in))
0.208 (8.2)
0.663 (26.1)
Infinity

Poisson’s Mass density
3
ratio
(kg/m )
0.35
2300
0.40
2000
0.45
1800

Elastic Modulus
(MPa (ksi))
703.3 (102)
55.8 (8.1)
57.2 (8.3)

The modulus is more reasonable in trial II, although the deflection basin may generate a
higher RMSE. The deflection basins from the FWD test and the backcalculation results by
MODCOMP5 are shown in Fig. 28. It shows that the deflection basin is simulated reasonably
well.

58

Deflection Basin

Maximum Deflection (mils)

0
5
ModComp-1
FWD test
ModComp-2

10
15
20
25
30
35
0

10

20

30
40
Radius distance (in)
Sensor Location (in)

50

60

70

Figure 28: Deflection basin during the FWD test and two cases of MODCOMP5
simulation at LTPP station 06-0565 in 1999

Using laboratory measured values of 6 psi confinement, the subgrade resilient modulus is
83.2 MPa (12.1 ksi). The elastic modulus obtained by the backcalculation by MODCOMP5 is
57.2 MPa (8.3 ksi), and the modulus by the wave velocity method is 101.5 MPa (14.7 ksi). Thus,
the difference between the experimental data and the backcalculation is 31.3%, compared to
21.5% when using the proposed wave velocity method. These values are summarized in Table
17.
Table 17: Elastic modulus of the subgrade layer at station 06-0565 in 1999

Lab_ measured Predicted by
k-value
MR @ 6psi
(Ek-value)
confinement
83.2 MPa
101.5 MPa
(12.1 ksi)
(14.7 ksi)

MODCOMP5
backcalculation
(EMODCOMP)
57.2 MPa
(8.3 ksi)
59

Difference
between MR

Difference
between MR and

and Ek-value

EMODCOMP

21.5 %

31.3 %

3.5.4 Discussion of Case Studies Obtained from LTPP
The three examples illustrate that the proposed method generates reasonable results. The
results obtained by wave propagation method are closer (20% to 80%) to the experimental data
than the backcalculation results obtained by MODCOMP5 (30% to 300%). However, this
observation is based on the lab-measured resilient modulus at the confinement of 6 psi, which
may not be right for the case studies above. Another experimental model is applied to estimate
the field elastic modulus of each subgrade layer, based on the possible compaction level during
construction and the load applied by FWD test.

3.5.4.1 Effect of Confinement
From experimental data, the resilient modulus (kPa), or elastic modulus, of the subgrade
layer is usually expressed as the Uzan model (Chen et al., 2007), in Eq. (29):

M R = k1θ k2τ k3 ,
where: θ is the bulk stress (kPa);

τ=

1
3

(29)

(σ 1 − σ 2 )2 + (σ 1 − σ 3 )2 + (σ 3 − σ 2 )2

is the octahedral shear stress (kPa), and

σ1, σ2 ,

and

σ3 are the principal normal

stresses; and k1, k2, k3 are the regression coefficients.
If the axial stress is almost identical to the confinement pressure, the octahedral shear stress
is around zero, and the regression is not easy, so the modified model, Cornell model (Irwin,
1994), is also implemented to better capture the above relationship, as in Eq. (30).

60

M R = k1θ

k2

(τ + 1)

k3

(30)

The experimental results of the above three examples are summarized for different cases of
axial load and confinement based on the LTPP database, and the RMSE technique is applied to
find the regression coefficients. There are no laboratory data for station 06-0565; however, the
experimental data for 06-0564 are available with two locations (BA9, BA10). One can assume
the subgrade should be the same since the two stations are separated for only 100 ft. The results
of data regression are shown in Table 18.

Table 18: Results of data regression of resilient modulus on stress parameters using Uzan
and Cornell Models

Station ID
04-1036
32-0101
06-0564-BA9
06-0564-BA10

Uzan model

k1

29295
4991
45994
24216

k2

0.1963
0.5341
0.1444
0.2177

k3

-0.0278
-0.0708
-0.0322
-0.0493

2

R

0.631
0.597
0.364
0.277

Cornell model

k1

29696
5224
44761
26140

k2

0.2068
0.5427
0.1690
0.2240

k3

-0.0521
-0.1027
-0.0665
-0.0910

R2
0.712
0.614
0.424
0.331

With the exception of data from station 06-0564 at location BA10, all other data are
2

reasonable with regression R >0.40 (Baladi, 2010). For simplification and better estimation, the
experimental data from 06-0564-BA9 were used for station 06-0565 for further analysis.
The actual confinement of the subgrade layer is unknown, although it was assumed to be
41.3 kPa (6 psi) in the design procedure. The coefficient of lateral earth pressure at rest (K0) for
subgrade may vary from 0.4 to 2.0, depending on the compaction level during construction. The
overburden stress in the Z-direction can be approximately calculated based on the mass density
of each layer, as the thickness of each layer is known. The stress status inside the subgrade layer
61

can be calculated based on the coefficient of lateral earth pressure at rest (K0).
For the FWD test, the stress inside the pavement can be directly obtained by the layered
elastic program JULEA, and the total stress is the summation of the stress at rest added to the
stress from the FWD load. Since the subgrade layer is semi-infinite, only the stress at 1 ft below
the subgrade surface (the compaction depth in construction) is considered, and the horizontal
distances are chosen for the location of sensors 5, 6 and 7.
The calculation variables and outcomes at station 04-1036 are shown in Table 19, given the
maximum FWD test pressure is 749 kPa (108.63 psi).

Table 19: The stress status of the subgrade layer in the FWD test for different K0 values

Sensor Horiz.
No. dist. (in)
5
6
7

24.0
36.0
60.0

Depth At rest FWD FWD FWD
(in) σ (psi)σ (psi) σ (psi) σ (psi)
z
x
y
z
35.4 2.52
-0.14
0.22
1.54
35.4 2.52
-0.10
0.35
0.82
35.4 2.52
-0.04
0.27
0.20

Octa.
k0=0.4 k0=1.0 k0=2.0 τ (psi)
6.2
9.2
14.2
0.98
5.6
8.6
13.7
0.67
5.0
8.0
13.0
0.28
Total Bulk Stress (psi)

Based on the bulk stress and octahedral shear stress ranges, the range of the resilient
modulus (MR) can be determined, given the regression coefficients in Table 18. A similar
procedure was done for stations 32-0101 and 06-0565. The range of MR calculated using the
non-linear model, the elastic modulus from forward calculation by Rayleigh wave method, and
backcalculation by MODCOMP5 for the three cases in the LTPP database are shown in Table 20.

62

Table 20: The range of elastic moduli predicted by nonlinear models, measured in the lab,
estimated by the proposed method, and backcalculated by MODCOMP5 for the case
studies (MPa)

Lab-Measured
Estimated by the MODCOMP5
MR @ 6 psi proposed method Backcalculation
Minimum Average Maximum
confinement
04-1036
58.4
64.6
72.7
102.0
156.9
222.7
32-0101
43.9
59.8
79.0
64.0
120.6
281.0
06-0565
79.0
86.2
93.8
83.2
101.5
57.2
The difference in subgrade moduli predicted by the three methods was smallest at station
Station
ID

Nonlinear model prediction

06-0565, and largest at station 04-1036. As the bulk stress affects the MR more than the
octahedral shear stress, the bulk stress, or the confinement is more important in determining the
MR of the subgrade layer. The maximum bulk stress at station 04-1036 is 13 psi or less, as shown
in Table 16, meaning the confinement is much less than 6 psi. At station 06-0565, the bulk stress
is much higher than that at 04-1036 (values are not shown here), as the pavement structure is
thicker in the base layer, and the confinement pressure may be around 6 psi.
The confinement of the subgrade layer is unknown, and the error may be less if the actual
confinement is known. The Poisson’s ratio and mass density are assumed in the proposed method,
and that have a linear effect on the final result. In the above examples, the mass density is
estimated from Stubstad (2002). If the mass density were more accurate, the error would be
smaller.

3.5.4.2 Effect of Sampling Intervals
Theoretically, it is possible to accurately predict the modulus of subgrade layers, if the
sampling interval is small enough and more sensors are available. Practically, the sampling time
is around 0.4 ms, and there are typically 7 sensors within 72 inches from the FWD load, so the
error may be significant. Assuming the S-wave velocity is 150 m/s, and R-wave velocity is 142.4
63

m/s, the time required for each wave to arrive at each sensor in the FWD test is listed in Table
21.

Table 21: Arrival time for S-wave and R-wave for each sensor due to the FWD test

Sensor
Location (m)
S-wave (ms)
R-wave (ms)

1
0
0
0

2
0.3
2.0
2.1

3
0.6
4.0
4.2

4
0.9
6.0
6.3

5
1.2
8.0
8.4

6
1.5
10.0
10.5

7
1.8
12.0
12.6

If the sampling interval is 0.4 ms in the FWD test, the time difference between S-wave and
R-wave propagation cannot be clearly identified for close sensors, which causes a significant
error for subsequent sensors. The vibrations of R-wave and S-wave are in opposite directions
(Fig. 9 in section 3.2.1). Because the two waves interact, the occurrence time of the peak
deflection is not accurate, and thus the error can be as high as 30%. The error difference between
numerical examples and case studies from the LTPP can be explained by differences in sample
time intervals. In order to reduce the error in estimation of subgrade elastic modulus from wave
propagation, the frequency of sampling should be increased.

3.6 Numerical Example of Bedrock under Subgrade layer
As shown in section 3.2.4, k-values vary if the elastic modulus of the subgrade layer is
higher than 500 MPa (70 ksi). However, in some cases, a shallow bedrock may underlie the
subgrade. In this section, the effect of existing shallow bedrock is investigated by numerical
simulation.
Bedrock is assumed to be located 5 ft below the surface of the subgrade layer, and the
time-history is simulated by SAPSI. Table 22 is a summary of the material properties of the

64

pavement profile.
Table 22: Material properties of a three-layer pavement structure with bedrock

Physical
Layer
AC
Base
Subgrade
Bedrock

Elastic Modulus
(MPa (ksi))
Experimental data
108 (21.8)
34.5 (5.0)
7000 (1000)

Poisson’s
ratio
0.35
0.30
0.45
0.25

Mass density
3
(kg/m )
2300
2000
1500
2800

Thickness
(m (in))
0.15 (6)
0.3 (12)
1.5 (60)
Infinity

The FWD load is assumed as a haversine function with a duration time of 35 ms, and the
maximum pressure is 707 kPa with a loading plate radius of 0.15 m. The simulated FWD
time-history is plotted in Fig. 29, and the time delay of each sensor is given in Table 23.

Deflection (m)

Deflection history with Bedrock (5 ft)
7.E-04
6.E-04
5.E-04
4.E-04
3.E-04
2.E-04
1.E-04
0.E+00
-1.E-04
-2.E-04

sensor 1
sensor 3

0.04

sensor 6

sensor 7

0.02

sensor 4

sensor 5

0

sensor 2

sensor 8

0.06

0.08

0.1

time (sec)
Figure 29: SAPSI simulated time-history of FWD test with bedrock underneath

65

Table 23: Time delay of each sensor by SPASI simulation with bedrock

Sensor
1
Distance (m)
0.00
Time delay (s) 0.0053

2
0.20
0.0056

3
0.30
0.0059

4
0.46
0.0064

5
0.61
0.0069

6
0.91
0.0083

7
1.22
0.01

8
1.52
0.0122

Fig.28 shows there is free vibration after the FWD test, due to the stiff layer underneath. The
time delay vs. location is plotted in Fig. 30.

Time delay vs. location (Bedrock 5 ft)

Time Delay (s)

0.014
0.012
0.01
0.008
0.006

y = 5.7743E-03x + 3.1900E-03

0.004

R2 = 9.8965E-01

0.002
0
0.00

0.40

0.80

1.20

1.60

Location (m)

Figure 30: Time delay vs. location for a numerical example of FWD test with bedrock

The slope can be obtained by linear regression, as m = 0.005774 s/m. Following the
proposed procedure in section 3.3, the elastic modulus of the subgrade layer can be calculated as:
E = 368.0 MPa, which is quite different from the modulus of the subgrade layer or the bedrock
layer. It is the combination of the two layers, as the more distant sensors are significantly

66

affected by the interaction at the interface of the subgrade layer and bedrock.
Therefore, the proposed method is not appropriate in the case of shallow bedrock. In wave
propagation theory, the path of wave propagation changes if there is a stiff layer underneath and
the wave is called Love wave instead of Rayleigh wave. A detailed explanation can be found in
reference (Richart et al., 1970). However, the existence of shallow bedrock can be easily
identified by the free vibration response of the pavement under FWD testing.

3.7 Summary
Both numerical examples and case studies obtained from the LTPP database showed that
the subgrade elastic modulus calculated by the proposed method based on wave
propagation theory was acceptable for pavement engineering application, compared to
the static backcalculation, no matter whether there is GWT underneath or not.
The estimation of the subgrade elastic modulus can be improved in field FWD tests by
increasing the sampling frequency of each sensor.
When shallow bedrock presents, the subgrade elastic modulus cannot be estimated
accurately, as the Rayleigh wave is theoretically based on one layer of semi-infinite
medium, and it is only valid for layered elastic systems with higher modulus of the
surface layers.
The existence of shallow bedrock underneath the subgrade layer is indicated by free
vibration in the pavement response due to the FWD test; in this case, the elastic modulus
of subgrade layer cannot be obtained by the proposed method.

67

Chapter 4. Layered Viscoelastic Forward
Solution
4.1 Introduction
In each loop of any backcalculation algorithm, there is a forward solution that predicts the
theoretical response, which is compared to the field measured response. Since a backcalculation
program often converges after many loops of iteration, time-efficiency and accuracy are critical
factors.
The forward viscoelastic solution is presented in this chapter. It utilizes the concept of
‘quasi-elastic’ approximation suggested by Schapery (1974). Both time-efficiency and accuracy
are compared, and the viscoelastic solution is chosen for further analysis, mainly based on the
consideration of time-efficiency. The error and limitation of the proposed solution is investigated
later in this chapter.

4.2 Layered Viscoelastic Forward Solution Algorithm
4.2.1 Theoretical Background
The time-dependent response of linear viscoelastic material subjected to a random loading
history can be computed using the following Boltzmann’s superposition integral (Schapery,
1974):

68

t

ve
RH (t − τ ) dI (τ )
∫

ve

R (t ) =

τ =0

where

R ve (t )

is the linear viscoelastic response,

unit step function of input ( I (t ) =

(31)
,

ve
RH (t )

is the viscoelastic response to a

H (t ) , where H(t) is the Heaviside step function) and

dI (τ )

is the change in input at time τ .
In the case of a layered system with a circular load as shown in Fig. 31, the time dependent
input is the contact stress (i.e., I(t) = σ(t)) and the response of interest may be the deflections (i.e.,

R ve (t ) = u ve (t ) ) at certain locations. Then we can rewrite Eq. (31) as follows:

t

u ve (t ) =

ve
u H (t − τ ) dσ (τ )
∫

τ =0
where

ve
uH

(32a)
,

is the deflection due to a unit contact stress (i.e., σ(t)=1). It is noted that

ve
uH

can be either the vertical or radial deflection, although in this paper, only the vertical surface
deflection is considered. The generalized formula for Eq. (32a) can be written, in cylindrical
coordinates, as follows:

t

ve
uvertical (t,z,r) =

∫u

τ =0

ve
H − vertical

(t − τ,z,r) dσ (τ )

(32b)

t

ve
uradial (t,z,r) =

∫u

τ =0

ve
H − radial

69

(t − τ,z,r) dσ (τ )

(32c)

ve
uvertical (t,z,r)

and

ve
uradial (t,z,r) , respectively, are the displacements in vertical

and radial directions, observed at time t and at location (r,z);

u ve−vertical (t − τ,z,r)
H

and

uve−radial (t − τ,z,r) , respectively, are the viscoelastic deflections due to a unit contact
H
stress (i.e., σ(t)=1).
r

σ

H
1
H2

H=Inf

Ε(t), ν1
Ε ,ν

AC
Base

2

Subgrade

2

Εn, νn

Figure 31: Typical geometry of a pavement structure
ve

The viscoelastic deflection due to a unit contact stress ( u H ) can efficiently and accurately
be computed by using Schapery’s ‘quasi-elastic’ approximation (Schapery 1965; 1974).
Quasi-elastic theory states that unit viscoelastic response (e.g.,

u ve (t) = surface deflection) at
H

any time t can be approximated by the unit elastic response (e.g.,

u e (E ve (t))
H

= elastic

surface deflection), which is calculated using the modulus ( E ) equal to viscoelastic modulus

E ve (t )

evaluated at time t, i.e.:

ve
u H (t) ≅ u e (E ve (t)) ,
H

70

(33)

where

e
uH

is the elastic deflection due to a unit step load computed using a layered elastic

solution where E

= E ve (t)

is utilized for the AC layer.

The detailed derivation of Eq. (33) can be found in Levenberg (2008) and will not be
repeated here for brevity. In this implementation, the unit response

ve
e
u H (t ) ≅ u H (t )

values

at the points of interest were computed using the CHEVLAY2 layered elastic analysis program.
Then the convolution integral in Eq. (32b) is used to calculate the viscoelastic deflection

u ve (t ) . The algorithm is described in the following section.
4.2.2 Algorithm Steps
Based on the theoretical formulation, the following steps are applied to implement the
Viscoelastic (VE) solution from the linear elastic layered solution of CHEVLAY2.
1. Define the geometry (layer thicknesses, contact pressure…etc) of a layered system
similar to the one in Fig. 31.
2. Select a stress versus time history, σ(t), and divide the data into Ns discrete intervals as
shown in Fig. 32.

71

120
Figure 2

Stress (psi)

100
80
60

Number of time intervals: Ns =50

40
20
0

0

0.02

0.04

0.06

0.08

0.1

t(s)

Figure 32: Discretization of stress history

3. Divide the relaxation modulus E(t) mastercurve into NE number of time steps (Fig. 33).

8

10

Figure 3
6

E(t) psi

10

4

10

Number of time intervals: NE=100
2

10
-6
10

-4

10

-2

10
t(s)

0

10

Figure 33: Discretization of the relaxation modulus mastercurve

72

2

10

4. Calculate the elastic response of the structure to a unit stress (step function of the load
input) using the E(ti) evaluated at different times (i.e., t1, t2, t3

….tNE). In this present

implementation, the surface deflections at several radial distances to a circular plate load
shown in Fig. 31 were of interest. Therefore, these surface deflections were computed
using the CHEVLAY2 program using the modulus value corresponding to different times
in Fig. 34; i.e., E(t1),

E(t2), E(t3), E(t4)… E(tNE):

ve
e
u H (ti ) ≅ u H

calculated using E(ti) where i=1,2,3 …NE

(34)

-1

10

rc = 0 in
uVE(t) = uE (t) (inches)
H
H

rc = 13 in
rc = 21 in

-2

10

rc = 35 in
rc = 49 in
rc = 63 in

-3

10

-4

10

-6

10

-4

-2

10

10
t(s)

0

10

2

10

Figure 34: Deflections calculated for points at different distances from the centerline of
the circular load at the surface

73

Fig. 34 shows the

e
uH

values calculated for points at different distances from the

centerline of the circular load at the surface. These curves are herein called “unit response
mastercurves”.
5. Calculate the viscoelastic response using the discrete form of Eq. (33) given in Eq. (35)
below. Eq. (35) is evaluated at each discrete time

Fig. 32. Fig. 35 below illustrates the

dσ (τ j )

ti

using the stress history shown in

in Eq. (35) for each time step

τj.

i

ve
u (t i ) = ∑ u H (t i − τ j ) dσ(τ j ),
ve

(35)

j= 0

i = 1, 2,...Ns

where

120

Stress (psi)

100
80
60
dσ (τ j )

40

δσ(τ)

20
0

τj

0

0.02

τ
Figure 35: Illustration of the

0.04

0.06

0.08

0.1

t(s)

dσ (τ j )

74

in Eq. (32b) for each time step

τj

In order to illustrate an example, the viscoelastic surface deflections of the three-layered
pavement structure shown in Fig. 36 are computed. Fig. 37 shows the vertical surface deflections
at points located at different radial distances from the centerline of the load, and clearly shows
the relaxation behavior of deflection at each point.
σ=102.6 psi

H1=3.94”
H2=11.8”

1

r=5.9”
2
3

6

5

E(t )= Fig. 15, ν1=0.35

AC
Base

H3=Inf Subgrade

4

E2=1567.5 psi ν2=0.35
E3=1451.4 psi, ν3=0.35

Figure 36: Example of a pavement structure used for VE forward calculation

0.15
rp1= 0 in
rp2= 13 in
rp3= 21 in

0.1
uVE(t) in

rp4= 35 in
rp5= 49 in
rp6= 63 in

0.05

0
0

0.02

0.04

0.06

0.08

0.1

t(s)

Figure 37: Examples of computed VE surface deflections at different radial distances from
the centerline of the load

75

4.3 Verification of the Layered Viscoelastic Forward Solution
First, verification of ‘quasi-elastic’ approximation was performed. Chen (2009) introduced a
theoretical solution for response of viscoelastic (VE) pavement layers. The unit responses
computed using Chen’s method (Chen and Pan, 2007) was compared against the ‘quasi-elastic’
approximation shown in Eq. (33). Second, the proposed VE solution is compared against the
well-known dynamic solutions by SAPSI and LAMDA. The effect of wave propagation was
eliminated from these dynamic solutions. Later, the time-efficiencies of both solutions are
compared in the next section.

4.3.1 Verification of Quasi-Elastic Approximation
A semi-analytical solution for a multilayered viscoelastic pavement under surface loading
was derived by Chen (2009; Chen et al, 2009). This solution was used to verify the accuracy of
the quasi-elastic approximation. The structural properties of the pavement used in this
verification are listed in Table 24.
Table 24: Material properties of the pavement structure used in the verification example

Physical Layer
AC
Base
Subgrade

Elastic Modulus
(MPa (ksi))
Experimental data
108 (15.6)
34.5 (5.0)

Poisson’s Ratio
0.35
0.30
0.45

Thickness
(m (in))
0.15 (6.0)
0.30 (12.0)
Infinity

Only 11 terms of the Prony series of E(t) can be implemented by the solution developed by
Chen (2009). The input E(t) curve was therefore approximated with a 11-term Prony series, as
shown in Fig. 38. It is noted that 11-term Prony series did not accurately approximate the actual
E(t) curve.

76

Implemented VS. actual E(t) of one AC layer

E(t) (MPa)

1.E+05

Actual

1.E+04

Implemented

1.E+03
1.E+02
1.E+01
1.E+00
1.E-15

1.E-10

1.E-05

1.E+00

1.E+05

1.E+10

1.E+15

Reference time (s)
Figure 38: Implemented and actual E(t) for one AC mixture

Comparison between analytical and quasi-elastic solution for
pavement response under step load
4.E-02
deflection (m)

3.E-02
3.E-02
2.E-02
2.E-02
1.E-02
5.E-03
0.E+00
0.E+00

sensor 1
sensor 5
1.E-02

sensor 2
sensor 6
2.E-02

3.E-02

sensor 3
Sensor 7
4.E-02

5.E-02

sensor 4

6.E-02

7.E-02

time (ms)

Figure 39: Implemented and actual E(t) for one AC mixture
(Solid lines represent the analytical solution, symbols represent the quasi-elastic
approximation)

77

Fig. 39 illustrates that the quasi-elastic approximation did a very good job in predicting unit
response. The maximum error was around 3% for sensor 1, and error decreased as sensor number
increased. This difference might be due to the 11-term Prony series approximation of E(t) used in
Chen’s method.

4.3.2 Comparison between the Layered VE Solution and Dynamic Solutions
The layered VE forward solution is also verified against the well-known dynamic solutions,
SAPSI and LAMDA. As discussed previously, wave propagation is ignored in the layered VE
solution. Therefore, the deflection time histories given by SAPSI and LAMDA were shifted, so
that the peaks of all sensors coincide.
In order to illustrate the absence of wave propagation, a FWD load pulse of duration 35 ms
was hypothetically apply to an infinitely deep AC layer, and the response was simulated using
the layered VE solution. The time-histories of the pavement response due to the FWD test
simulation are shown in Fig. 40.

78

FWD deflection history

Deflection (in)

1.2E-03
0 in
12 in
24 in
48 in
72 in

1.0E-03
8.0E-04
6.0E-04

8 in
18 in
36 in
60 in

4.0E-04
2.0E-04
0.0E+00
0.00

0.01

0.02

0.03

0.04

0.05

0.06

0.07

Time (s)
Figure 40: FWD test simulation by layered VE solution for a half-space AC layer

Delayed recovery of the displacements/deflection, which is a clear indication of viscoelastic
response, is visible in Fig. 40. However, there is no time delay of the peaks between the sensors,
i.e., the wave propagation cannot be simulated by layered VE solution. Since the focus of this
study is on the viscoelasticity property, not the dynamic response of the system, the wave
propagation response is eliminated by shifting the time histories of all sensors such that their
peaks occur at the same time. .
The layered VE solution is superimposed on the time-shifted dynamic solutions, in order to
better compare the shape of the deflection time histories, and to show the accuracy of the layered
VE solution.
Table 25 shows that the properties of a pavement structure used in comparing dynamic and
79

layered VE solutions. Fig. 41 and Fig. 42 show the response of the pavement under a FWD test
calculated by SAPSI and LAMDA, respectively.
Table 25: Properties of a pavement used in comparison of dynamic and VE solutions

Physical Layer
AC
Base
Subgrade

Elastic Modulus
(MPa (ksi))
Experimental data
108 (15.6)
34.5 (5.0)

Poisson’s Ratio
0.35
0.40
0.45

Density
3
(kg/m )
2300
2000
1500

Thickness
(m (in))
0.10 (4.0)
0.3 (12.0)
Infinity

FWD deflection history
0.700

r=0 m
r=0.33
r=0.52
r=0.88
r=1.25
r=1.61

Deflection (mm)

0.600
0.500
0.400
0.300
0.200

m
m
m
m
m

0.100
0.000
-0.100
0.000

0.020

0.040

0.060

0.080

0.100

Time (second)
Figure 41: Sensor deflection time histories predicted by SAPSI

80

0.120

FWD deflection history
r=0 m
r=0.33
r=0.52
r=0.88
r=1.25
r=1.61

0.700
Deflection (mm)

0.600
0.500
0.400
0.300
0.200

m
m
m
m
m

0.100
0.000
-0.100
0.000

0.020

0.040

0.060

0.080

0.100

0.120

Time (second)
Figure 42: Sensor deflection time histories predicted by LAMDA

The layered VE solution and the time-shifted dynamic solutions for the same example are
shown in Fig. 43 and Fig. 44. As shown in these figures, layered VE solution and dynamic
solutions by SAPSI and LAMDA match very well. It is noted that the computational efficiency
of VE solution is much better than SAPSI and LAMDA, where it takes around 20 minutes to
obtain a solution from SAPSI and LAMDA, while it only takes 1 minute, or even a few seconds,
to reach a solution using layered VE solution.

81

Deflection (mm)

VE solution VS. SAPSI

VE Sensor 0
VE Sensor 2
VE Sensor 4
SAPSI Sensor 0
SAPSI Sensor 2
SAPSI Sensor 4

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
0

0.01

0.02

0.03

0.04
0.05 0.06
Time (second)

0.07

VE Sensor 1
VE Sensor 3
VE Sensor 5
SAPSI Sensor 1
SAPSI Sensor 3
SAPSI Sensor 5

0.08

0.09

0.1

Figure 43: Comparison of deflection time histories from VE and SAPSI solutions
(Solid lines represent the VE solution, symbols represent the SAPSI result)

VE solution VS. Lamda

Deflection (mm)

VE sensor 0
VE sensor 2
VE sensor 4

0.02

0.04

Lamda Sensor 1
Lamda Sensor 3

Lamda Sensor 4

0

VE sensor 3
VE sensor 5

Lamda sensor 0
Lamda Sensor 2

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1

VE sensor 1

Lamda Sensor 5

0.06

0.08

Time (second)
Figure 44: Comparison of deflection time histories from VE and LAMDA solutions
(Solid lines represent the VE solution, symbols represent the LAMDA result)

82

0.1

4.3.3 Verification Examples using Pavement Structures from the LTPP
Database

Table 26: Several examples of pavement structure in the SPS-1 project

Case No. Layer No.
4
3
116
2
1

Layer Description
Original Surface Layer (Layer Type:AC)1.8 Inch
AC Layer Below Surface (Binder Course) (Layer Type:AC)2.1 Inch
Base Layer (Layer Type:TB)12 Inch
Subgrade (Layer Type:SS)

Layer No.
6
5
4
3
2
1

Layer Description
Original Surface Layer (Layer Type:AC)1.7 Inch
AC Layer Below Surface (Binder Course) (Layer Type:AC)1.5 Inch
AC Layer Below Surface (Binder Course) (Layer Type:AC)3.2 Inch
Base Layer (Layer Type:TB)5.2 Inch
Base Layer (Layer Type:GB)4 Inch
Subgrade (Layer Type:SS)

Layer No.
5
4
3
2
1

Layer Description

117

120

Original Surface Layer (Layer Type:AC)1.8 Inch
AC Layer Below Surface (Binder Course) (Layer Type:AC)1.8 Inch
Base Layer (Layer Type:PATB)4 Inch
Base Layer (Layer Type:GB)8 Inch
Subgrade (Layer Type:SS)

Layer No. Layer Description

123

6
5
4
3
2
1

Original Surface Layer (Layer Type:AC)1.8 Inch
AC Layer Below Surface (Binder Course) (Layer Type:AC)2 Inch
AC Layer Below Surface (Binder Course) (Layer Type:AC)2.4 Inch
Base Layer (Layer Type:TB)8 Inch
Base Layer (Layer Type:PATB)4 Inch
Subgrade (Layer Type:SS)

83

In order to further verify the layered VE solution, four pavement structures are selected from
the LTPP database. Table 26 shows the selected pavement structures from the SPS-1 experiment
of the LTPP database.
As shown in Table 26, there are several sublayers of AC. For computational efficiency,
similar layers were combined into one layer, as shown in Table 27.

Table 27: Layer information for each pavement structure

Case
No.
116

Physical
Layer
AC
Base
Subgrade

Elastic Modulus
Thickness
(MPa (ksi))
(m (in))
Experimental Data 0.099 (3.9)
200 (29.0)
0.305 (12.0)
100 (14.5)
Infinity

117

AC
Experimental Data
Base
200 (29.0)
Granular Base
150 (21.8)
Subgrade
100 (14.5)

0.163 (6.4)
0.132 (5.2)
0.102 (4.0)
Infinity

120

AC
Experimental Data 0.091m (3.6)
PATB*
180 (26.1)
0.102 (4.0)
Granular Base
150 (21.8)
0.204 (8.0)
Subgrade
100 (14.5)
Infinity

AC
Experimental Data 0.157 (6.2)
Base
200 (29.0)
0.204 (8.0)
123
PATB*
180 (26.1)
0.102 (4.0)
Subgrade
100 (14.5)
Infinity
Note: PATB = Permeable Asphalt Treated Base

Poisson’s
Ratio
0.35
0.40
0.45

Mass density
(kg/m3)
2300
2000
1800

0.35
0.40
0.40
0.45

2300
2000
2000
1800

0.35
0.40
0.40
0.45

2300
2000
2000
1800

0.35
0.40
0.40
0.45

2300
2000
2000
1800

The deflection time histories predicted by SAPSI, LAMDA (after eliminating the time delay,
or time-shifted dynamic solution) and by the VE solution are shown in Figures 45 to 48.

84

VE solution VS. SAPSI for Case 116
VE Sensor 0

VE Sensor 1

0.5
Deflection (mm)

0.6

VE Sensor 2

VE Sensor 3

VE Sensor 4

VE Sensor 5

SAPSI Sensor 0

SAPSI Sensor 1

0.3

SAPSI Sensor 2

SAPSI Sensor 3

0.2

SAPSI Sensor 4

SAPSI Sensor 5

0.4

0.1
0
-0.1
0

0.02

0.04
0.06
Time (second)

0.08

0.1

(a) VE solution and SAPSI
VE solution VS. Lamda for Case 116
VE Sensor 0

VE Sensor 1

VE Sensor 2

VE Sensor 3

VE Sensor 4

VE Sensor 5

0.4

Lamda Sensor 0

Lamda Sensor 1

0.3

Lamda Sensor 2

Lamda Sensor 3

0.2

Lamda Sensor 4

Lamda Sensor 5

0.6

Deflection (mm)

0.5

0.1
0
-0.1
0

0.02

0.04
0.06
Time (second)

0.08

0.1

(b) VE solution and LAMDA
Figure 45: Time-shifted dynamic solution and VE solution for case 116
(Solid lines represent the VE solution, symbols represent the SAPSI / LAMDA result)

85

VE solution VS. SAPSI for Case 117
VE Sensor 0

VE Sensor 1

0.4

VE Sensor 2

VE Sensor 3

0.35

VE Sensor 4

VE Sensor 5

0.3

SAPSI Sensor 0

SAPSI Sensor 1

0.25
0.2

SAPSI Sensor 2

SAPSI Sensor 3

SAPSI Sensor 4

SAPSI Sensor 5

Deflection (mm)

0.45

0.15
0.1
0.05
0
-0.05
0

0.02

0.04
0.06
Time (second)

0.08

0.1

(a) VE solution and SAPSI
VE solution VS. Lamda for Case 117
VE Sensor 0

Deflection (mm)

VE Sensor 3

VE Sensor 4

VE Sensor 5

Lamda Sensor 0

Lamda Sensor 1

Lamda Sensor 2

Lamda Sensor 3

Lamda Sensor 4

0.35
0.3
0.25
0.2
0.15

VE Sensor 1

VE Sensor 2

0.45
0.4

Lamda Sensor 5

0.1
0.05
0
-0.05
0

0.02

0.04
0.06
Time (second)

0.08

0.1

(b) VE solution and LAMDA
Figure 46: Time-shifted dynamic solution and VE solution for case 117
(Solid lines represent the VE solution, symbols represent the SAPSI / LAMDA result)

86

VE solution VS. SAPSI for Case 120
VE Sensor 0

Deflection (mm)

0.2

SAPSI Sensor 3

SAPSI Sensor 4

0.3

SAPSI Sensor 1

SAPSI Sensor 2

0.4

VE Sensor 5

SAPSI Sensor 0

0.5

VE Sensor 3

VE Sensor 4

0.6

VE Sensor 1

VE Sensor 2

0.7

SAPSI Sensor 5

0.1
0
-0.1
0

0.02

0.04
0.06
Time (second)

0.08

0.1

(a) VE solution and SAPSI
VE solution VS. Lamda for Case 120
VE Sensor 0

VE Sensor 1

0.5
Deflection (mm)

0.6

VE Sensor 2

VE Sensor 3

0.4

VE Sensor 4

VE Sensor 5

Lamda Sensor 0

Lamda Sensor 1

Lamda Sensor 2

Lamda Sensor 3

Lamda Sensor 4

Lamda Sensor 5

0.3
0.2
0.1
0
-0.1
0

0.02

0.04
0.06
Time (second)

0.08

0.1

(b) VE solution and LAMDA
Figure 47: Time-shifted dynamic solution and VE solution for case 120
(Solid lines represent the VE solution, symbols represent the SAPSI / LAMDA result)

87

VE solution VS. SAPSI for Case 123
VE Sensor 0

VE Sensor 1

VE Sensor 2

VE Sensor 3

VE Sensor 4

VE Sensor 5

0.3

SAPSI Sensor 0

SAPSI Sensor 1

0.25
0.2

SAPSI Sensor 2

SAPSI Sensor 3

SAPSI Sensor 4

SAPSI Sensor 5

0.45
0.4
Deflection (mm)

0.35

0.15
0.1
0.05
0
-0.05
0

0.02

0.04
0.06
Time (second)

0.08

0.1

(a) VE solution and SAPSI
VE solution VS. Lamda for Case 123
VE Sensor 0

VE Sensor 1

0.4
0.35
Deflection (mm)

0.45

VE Sensor 2

VE Sensor 3

0.3
0.25
0.2

VE Sensor 4

VE Sensor 5

Lamda Sensor 0

Lamda Sensor 1

Lamda Sensor 2

Lamda Sensor 3

Lamda Sensor 4

Lamda Sensor 5

0.15
0.1
0.05
0
-0.05
0

0.02

0.04
0.06
Time (second)

0.08

0.1

(b) VE solution and LAMDA
Figure 48: Time-shifted dynamic solution and VE solution for case 123
(Solid lines represent the VE solution, symbols represent the SAPSI / LAMDA result)

88

Figures 45 - 48 show that the typical difference between the results by the VE solution and
those by time-shifted dynamic analysis is about 6%, with the maximum difference as high as
10%. The difference appears to be slightly higher when the AC layer is thicker, or when there are
more physical layers in the pavement structure. Most importantly, the shapes of the response
pulses from dynamic and VE solutions are consistent with each other. Therefore, the four
numerical examples show that the difference between the responses from the dynamic solutions
(SAPSI or LAMDA) and from the VE solution is minimal.

4.3.4 A Numerical Example with the Bedrock
The existence of bedrock may cause problems in numerical simulation. The previous chapter
shows that the proposed method to estimate the elastic modulus of the subgrade layer by
Rayleigh wave is invalid if there is shallow bedrock underneath. In addition, the reflection of the
waves and resulting response can overshadow the effects of viscoelasticity. Fig. 49 shows the
responses calculated by SAPSI and layered VE solution. As expected, dynamic response by
SAPSI is significant, while layered VE solution does not show the effect of wave propagation.

89

Comparison of deflection history with
Bedrock (5 ft)

Deflection (m)

8.E-04

Sensor
Sensor
Sensor
Sensor

6.E-04
4.E-04

1
3
5
7

Sensor
Sensor
Sensor
Sensor

2
4
6
8

2.E-04
0.E+00
-2.E-04
0

0.02

0.04

0.06

0.08

0.1

time (s)
Figure 49: Time-shifted SAPSI and VE solution for shallow bedrock
(Solid lines represent the SAPSI result, symbols represent the VE solution)

Fig. 49 illustrates that the layered VE solution cannot accurately simulate the pavement
response if there is shallow bedrock underneath. Therefore, the layered VE solution should not
be used for simulating the pavement response if there is a shallow, stiff layer underneath. The
magnitude of error of the simulation depends on the relative level of stiffness between physical
layers, and the relative depth of the bedrock
A sensitivity analysis for the depth of the bedrock is done for the numerical example, and the
only input variable is the thickness of the subgrade layer, or the depth of the bedrock, and the
output variable is the relative error of the maximum deflection for each sensor between dynamic
solution SAPSI and VE forward solution. The result is shown in Fig. 50. Depending on the

90

accuracy of the simulation requirement, the effect of bedrock is insignificant if the depth is 15 ft
or more.

Relative Error of
Maximum Deflection

Relative Error vs. Bedrock Depth
70%
60%
50%
40%
30%
20%
10%
0%
-10%

Sensor
Sensor
Sensor
Sensor

0

5

10

1
3
5
7

15

Sensor
Sensor
Sensor
Sensor

2
4
6
8

20

Bedrock Depth (ft)
Figure 50: Sensitivity analysis for relative error of maximum deflection between VE
solution and SAPSI with bedrock depth

4.4 Time Efficiency for Layered VE Forward Solution
One of the primary reasons for implementing Schapery’s “quasi-elastic” approximation is its
extreme computational efficiency. It takes more than 25 minutes to calculate the pavement
response under step load by semi-analytical solution, and it takes around 20 minutes to simulate
the pavement response under FWD test by dynamic solution SAPSI or LAMDA, while it takes
less than one minute to calculate the response by VE solution. Many repetitions are required for
the forward simulation, so time efficiency is a very important factor.
Using a Pentium (4) 3.20 GHz computer with 1.99 GB ram, the computation time of the
results shown in Fig. 37 is 44.72 seconds, where the solution was for the 3-layered system shown

91

in Fig. 31 and NS = 50, NE=50. Table 28 shows the computation times for different numbers of
discrete time steps for the 3-layered system shown in Fig. 35. The difference in Table 28 is
defined as square root of the relative difference of the peak deflection compared to the most
difficult (reference) case with NS = 200 and NE = 200. Mathematically the difference is
expressed as:
2

1 N ⎛ ωi ( peak in case ) − ωi ( peak reference ) ⎞
⎟ ,
∑⎜
⎜
⎟
ωi ( peak reference )
N i =1 ⎝
⎠

Difference =

(36)

where: N = No. of Sensors.

Table 28: The computation times for different numbers of discrete time steps

Ns

NE

Elapsed time for computation (seconds)

Difference (%)

50
24
50
100
100
200

50
100
100
100
200
200

22.60
45.41
44.72
45.35
89.38
89.73

0.069
0.360
0.077
0.011
0.018
0.000

Table 28 shows that NS = 100 and NE = 100 is enough to generate a reasonable result, i.e.,
the sampling time for the FWD test is around 0.5 ms, and the time cost of each forward
simulation is around 45 seconds.
Therefore, the VE forward solution is accurate and time efficient, and is chosen as the
numerical simulation of the pavement under FWD test in future analysis.

92

4.5 Summary
A computation procedure for layered VE solution has been provided in this chapter. It
was simple in understanding, efficiency in calculation, and easy for implementation in
MATLAB.
The accuracy of the proposed forward solution was checked against a semi-analytical
solution for a step load, and was checked against dynamic solutions SAPSI and
LAMDA for FWD test simulation. The time delay of each sensor cannot be simulated,
and the deflection time-histories can be calculated within error around 5% for maximum
deflection.
The proposed layered VE forward solution was chosen for future backcalculation
because of its accuracy and time-efficiency.
If there is shallow bedrock underneath the pavement, the effect of mass inertia cannot be
ignored, and the proposed VE solution is not appropriate.

93

Chapter 5. Backcalculation Algorithm
5.1 Introduction
In this chapter, the outputs of the backcalculation program are identified, and the boundaries
of each output variable are measured from the lab as guidance for the convergence of the
program. Then the basic theory of inverse problem is briefly introduced. A detailed flowchart of
the backcalculation program is presented, and one numerical example of the backcalculation
result is shown in the end.

5.2 Variable Identification
The input variables are thickness of each layer, Poisson’s ratio of each layer, applied load
history and deflection time history of each sensor during a FWD test. The output variables are
the elastic modulus of base and subgrade layers, and a function of relaxation modulus E(t) or
dynamic modulus mastercurve |E*| of the AC layer. The elastic modulus of base and subgrade
layers is noted separately as Ebase and Esubgrade.
Mathematically E(t) is commonly expressed as Prony series with more than 30 variables,
however, it is almost impossible for any backcalculation to converge with so many variables. A
pre-described function must be required to reduce the unknown parameters. The sigmoid
function is preferred for |E*| in frequency domain and E(t) in time domain, as Eq. (1) for |E*|
and Eq. (2) for E(t) respectively. The backcalculation is used to calculate the parameters of the

94

sigmoidal function (a, b, c, d) for |E*| in frequency domain or E(t) in time domain respectively,
although the forward calculation is done in time domain.

5.3 Variable Boundaries
A reasonable range for the value of each output variable helps to converge the
backcalculation. If the intermediate output value crosses the boundary, it is forced to be within
the range by the algorithm.
The elastic modulus of base and subgrade layer is suggested to be less than 700 MPa (100
ksi). If not, it is not appropriate for the proposed wave propagation method to estimate the elastic
modulus of the subgrade layer (as discussed in Chapter 3) and for the proposed layered VE
solution to simulate the pavement response.
The boundaries of the parameters for the AC layer in the sigmoid function are found in the
experiment data of many mixtures. The seed value of the four parameters can be set as the
average value, and the boundary limit is identified via the minimum or maximum value of each
variable. The corresponding |E*| or E(t) components are discussed in the following separate
sections.

5.3.1 Boundary Limits for |E*|
About 30 sets of |E*| were obtained from FHWA (Kutay, 2008) and other sources (Mogawer
et al, 2010). Their coefficients for the dynamic modulus mastercurve |E*| are shown in Table 29.

95

Table 29: The |E*| parameters for different AC mixtures

Mixture Name
PG64-28 Control no-PPA (AI)
ALF Control 70-22
SBS 64-40
Terpolymer
CR-AZ
SBS LG
CR-TB
Air Blown
Fiber
Control 70-22 KIM
Advera
Control
Sasobit
PPA
PPA + Elvaloy
SBS
SBS + PPA
SBS PG64-34
PG64-28 Control no-PPA (AI)
PG64-28 with PPA (Hudson)
PG64-34 (SEM)
PG76-22 (Citgo)
PG64-22 with 12% GTR (Gorman)
PG64-28 no-PPA(AI)with2.0%Latex
WAM Control
WAM Foam
Asphamin
Control
Sasobit-1
Sasobit-2
Evatherm-1
Evatherm-2

a
1.699
0.797
1.338
1.468
-0.202
1.341
0.807
0.553
0.504
1.585
2.156
0.458
-0.26
1.552
1.579
1.72
1.568
0.851
1.569
1.49
1.689
1.646
1.245
1.507
0.92
1.13
1.261
1.195
1.423
1.266
1.112
1.191

b
2.765
3.593
2.843
2.799
4.695
2.857
3.519
3.868
3.838
2.955
1.801
3.772
4.907
2.408
2.317
2.21
2.505
2.391
2.65
2.677
2.394
2.476
2.851
2.678
3.371
3.168
2.985
3.127
2.791
3.012
3.162
3.144

c
0.823
1.372
0.164
0.491
1.122
0.876
1.066
1.285
1.535
1.077
0.546
0.948
0.787
0.602
0.582
0.519
0.464
0.157
0.855
0.777
0.156
0.954
1.151
0.905
1.387
1.244
1.284
1.333
1.423
1.375
1.245
1.258

d
0.483
0.503
0.561
0.583
0.346
0.556
0.407
0.405
0.504
0.558
0.854
0.389
0.268
0.702
0.543
0.73
0.643
0.682
0.668
0.639
0.681
0.705
0.483
0.572
0.595
0.641
0.598
0.552
0.578
0.543
0.546
0.497

The average, minimum and maximum values of each variable are listed in Table 30.

96

Table 30: The average and boundaries for the four parameters of |E*|

Average value
Minimum value
Maximum value

a
1.192
-0.26
2.156

b
3.016
1.801
4.907

c
0.93
0.156
1.535

d
0.563
0.268
0.854

5.3.2 Boundary Limits for E(t)
The relaxation modulus E(t) of the same 30 examples of AC mixture is expressed via the
sigmoid function, and values of the four parameters of the function are listed in Table 31.

97

Table 31: The E(t) parameters for different AC mixtures

Mixture Name
PG64-28 Control no-PPA (AI)
ALF Control 70-22
SBS 64-40
Terpolymer
CR-AZ
SBS LG
CR-TB
Air Blown
Fiber
Control 70-22 KIM
Advera
Control
Sasobit
PPA
PPA + Elvaloy
SBS
SBS + PPA
SBS PG64-34
PG64-28 Control no-PPA (AI)
PG64-28 with PPA (Hudson)
PG64-34 (SEM)
PG76-22 (Citgo)
PG64-22 with 12% GTR (Gorman)
PG64-28 no-PPA(AI)with2.0%Latex
WAM Control
WAM Foam
Asphamin
Control
Sasobit-1
Sasobit-2
Evatherm-1

a
1.721
0.841
1.342
1.472
0.045
1.351
0.895
0.682
0.561
1.598
2.155
0.558
0.150
1.553
1.586
1.720
1.570
0.851
1.571
1.493
1.689
1.647
1.278
1.515
0.934
1.136
1.272
1.216
1.438
1.289
1.132

b
2.735
3.540
2.832
2.790
4.400
2.842
3.411
3.717
3.772
2.937
1.800
3.645
4.378
2.405
2.306
2.208
2.501
2.388
2.646
2.671
2.391
2.473
2.810
2.666
3.352
3.159
2.971
3.101
2.772
2.984
3.137

c
0.326
0.860
-0.399
-0.095
0.731
0.313
0.634
0.855
1.025
0.512
-0.307
0.532
0.436
-0.100
0.030
-0.211
-0.180
-0.523
0.185
0.135
-0.523
0.247
0.653
0.326
0.789
0.601
0.681
0.774
0.837
0.823
0.692

d
-0.491
-0.515
-0.563
-0.584
-0.380
-0.559
-0.428
-0.431
-0.519
-0.562
-0.845
-0.411
-0.314
-0.699
-0.546
-0.725
-0.641
-0.679
-0.666
-0.638
-0.678
-0.702
-0.494
-0.574
-0.599
-0.642
-0.601
-0.559
-0.583
-0.551
-0.552

Evatherm-2

1.227

3.100

0.749

-0.508

The average, minimum and maximum values of each variable are listed in Table 32.

98

Table 32: The mean and boundaries for the four parameters of E(t)

Average value
Minimum value
Maximum value

a
1.234
0.045
2.155

b
2.964
1.800
4.400

c
0.356
-0.523
1.025

d
-0.570
-0.845
-0.314

5.4 Inverse Solution Theory
The objective of any backcalculation solution is to find a set of parameters. In this case, the
parameters are the set (a, b, c, d) of the AC mixture, the elastic modulus of base (Ebase), and the
elastic modulus of the subgrade layer (Esubgrade). In the end, the calculated deflection history
will match the measured values within a specified tolerance. To accomplish this goal, the
algorithm repeatedly adjusts the parameter values until a suitable match is obtained.

5.4.1 Basic Theory Background
The VE backcalculation algorithm in this research is an extension of the solution used in the
MICHBACK program (Harichandran et al., 1994). It uses the modified Newton method to obtain
a least squares solution of an over determined set of equations. In the MICHBACK solution,
these sets are real-valued and correspond to the peak deflection values. In this research, the
author will use deflection time histories, or many deflection basins (corresponding to different
times after eliminating the time delay of all sensors), since the proposed backcalculation scheme
uses a layered VE solution to predict the time-dependent deflection basins.
Ebase and Esubgrade are not known, and the AC layer can be expressed as a function of four
parameters (a, b, c, d). The unknown vectors become

99

{X } = {a

b c d

Ebase

Esubgrade }

T

, assuming all other information for base,

sub-base and roadbed layers are known.
The vector of measured responses is therefore expressed as:

⎧[w1 (t s ) ... wm (t s )],..., ⎫
{U (t )} = ⎪[w1 (t k ) ... wm (tk )],...,⎪ ,
⎬
⎨
⎪ w (t ) ... w (t ) ⎪
m f
⎭
⎩ 1 f
T

[

(37)

]

where m is the number of sensors in the FWD test, wj(ts) is the deflection of sensor j at the
starting time and wj(tf) is the deflection of sensor j at the final time of the specified range.
Following the derivation by Harichandran et al. (1994), the increment to the unknown
parameters in iteration i, {Δx}i , is obtained by solving the linear set of equations:

{ }

ˆ (t ) i + [G ]i {ΔX } = { (t )},
U
U
i

{Uˆ (t )}

(38)

i

where:

is the vector of deflections at individual time steps within the specified time

range, computed using the estimates of the parameters

ˆ
{x}i

at iteration i, and

gradient matrix at iteration i given by:

⎡ ∂ {U ( t )} ⎤
⎥
[G ] = ⎢
⎢ ∂ {X} ⎥ {X}={x}i
⎣
⎦
ˆ
i

100

(39)

[G ]i is the

The partial derivatives in the gradient matrix must be evaluated numerically using:

(

∂w j (t )
∂xk

ˆ
{X }={x}i

)

( )

ˆ}i − w j {x}i
ˆ
w j [R ]{x
=
,
i
ˆ
rx k

j = 1,2,..., m, k = 1,2,3,4,5,6

(40)
,

th

where: [R] is a diagonal matrix with the k diagonal element being (1+r) and all other diagonal
elements being one. A separate call to the forward calculation program is required to compute the
partial derivatives in each column of the gradient matrix.
If there are n individual time steps, Eq. (38) represents a set of m by n equations for 6
unknowns. Since there are more equations than unknown, a more robust method for solving the
problem is to use the singular value decomposition (SVD). This method will be briefly
introduced later.

After the increments

for

{X }

{Δx}i

are obtained by solving Eq. (38), the revised parameters

are obtained from Eq. (40), as below:

ˆ}i +1 = {x}i + {Δx}i
ˆ
{x

(41)

The iteration is terminated when the changes in the six parameters are smaller than a set of
specified tolerances:

101

⎛ xk+1 − xk
ˆi
ˆi
abs⎜
⎜ x i +1
⎝ ˆk

⎞
⎟≤ε
⎟
⎠

k = 1,2,3,4,5,6

(42)

5.4.2 SVD Method
SVD is a very powerful set of techniques for dealing with sets of equations or a matrix that
are either singular or numerically very close to singular. SVD methods are based on the
following theorem of linear algebra (Press et al, 1989). Any M×N matrix [A] whose number of
rows M is greater than or equal to its number of columns N, can be written as the product of an
M×N column-orthogonal matrix [U], an N×N diagonal matrix [W] with positive or zero elements
(the singular values), and the transpose of an N×N orthogonal matrix [V].
The matrix [A] can be decomposed as three matrices as following:

[A] = [U ] ∗ [W ] ∗ [V ]T

⎡ω1
⎢
ω2
⎢
= [U ]*
⎢
⎢
⎣

⎤
⎥
⎥ ∗ [V ]T
⎥
,
⎥
ωN ⎦

(43)

[U ]* [U ]T = I , [V ]* [V ]T = I , and [W] is diagonal matrix.
where:
For the following ill conditioned system of equations:

[A]* {X } = {b}
Inversing [A] by SVD method,

{X } can then be expressed as:
102

(44)

⎡1
⎢ω
⎢ 1
⎢
⎡1⎤
{X } = [V ]* ⎢ ⎥ * [U ]* {b} = [V ]* ⎢
⎣W ⎦
⎢
⎢
⎢
⎢
⎣
T
T
where: [U ] * [U ] = I , [V ] * [V ] = I , and

1

ω2

⎤
⎥
⎥
⎥
⎥ * [U ]* {b}
⎥
,
⎥
1 ⎥
ωN ⎥
⎦

(45)

⎡1⎤
⎢W ⎥ is a diagonal matrix.
⎣ ⎦

5.4.3 Truncating Singular Values
Reference (Press et al., 1989) defines the condition number of a matrix as a ratio of the
largest ωi of to the smallest ωj.

ωmax
Condition number =
ωmin

(46)

A matrix [A] is singular or ill conditioned if its condition number is too large. Reference
(Press et al., 1989) suggested this number should be adjusted according to the experiment with
the specific problem.
Define the threshold of truncation as:

ω j < ωmax *10−Threshold ,
j = 1,2,..., N
103

(47)

After selecting the threshold of a condition number, the SVD algorithm will simply

1
replace

ωj

with zeros.

To reduce the possibility of ill conditioning in the inverse problem, the technique of scaling
can be applied. In this program, the magnitude of Ebase and Esubgrade is within the magnitude of
100 MPa (14.5 ksi), but the magnitude of the four parameters in the sigmoid function for AC
layer is around 0 – 5, as shown in Table 30 & 32. The condition number would be very high, thus,
scaling technique is required. Ebase and Esubgrade are scaled by 100 MPa in this backcalculation
program. More information can be found in the references (Ji, 2005).

5.5 VE backcalculation of E(t)
Besides the above proposed Newton’s method, MATLAB internal function “fminsearch”
works very well for backcalculation problems, if the seed value is close to the actual solution.
Thus, the proposed Newton’s method is to scan all possible values and eliminate the local
minimum values for the error in Eq. (42).
Avoiding the local minimum problem in backcalculation is essential, so a random function
inside the backcalculation algorithm is used to prevent the solution converging to local minimum
th

values. A simplified flow chart is shown in Fig. 51, for the i loop of the backcalculation.

104

(i)

Input: X

Loop: k=1
k=k+1

(

)

⎡ abs X (i ) ( j ) − X (k ) ( j ) ⎤
(k )
Error = ∑ ⎢
⎥
X (i ) ( j )
⎢
⎥
⎣
⎦

2

Y
Error

N

(k)

<ε

k<i
N

Y
(i)

Random function: X

(i)

Output: X

Figure 51: Random Function to avoid local minimum problem

Fig. 51 shows the program cannot always avoid local minimum problem. However, as the
number of loops increases, the chance that the program will converge to the global minimum
point increases, at the cost of time-efficiency. The proposed wave propagation method in Chapter
3 is used to estimate the seed value of Esubgrade. A flow chart for this algorithm is shown in Fig.
52.

105

Input: Pavement structure profile,
FWD sensor profile, F(t), Wi(t)
Time delay ∆ti, deflection history Wi (tj)

Estimate subgrade Modulus by k-value

Seed value:

i=i+1

{X}

Boundary Check
(i)

Deflection basin Wi (tj) from VE solution
⎡
⎛ (i) ⎛ ⎞
⎛ ⎞⎞⎤
⎢ abs ⎜ w i ⎜ t j ⎟ − w i ⎜ t j ⎟ ⎟ ⎥
(i )
⎝ ⎠
⎝ ⎠⎠⎥
⎝
Error = ∑ ⎢
⎢
⎥
w ⎛t ⎞
i ⎜ j⎟
⎢
⎥
⎝ ⎠
⎣
⎦

(i)

Error <ε

2

Y

i<NMax

Avoid Local Minimum

Increment by Newton’s Method

N

Y
New Seed value:

{X}

MATLAB “fminsearch”

Output:

N

{X}

Figure 52: Algorithm of backcalculation procedure of E(t)

106

5.6 VE backcalculation of |E*|
An independent backcalculation algorithm is developed in this section if |E*| of the AC
mixture is concerned. The output variables are the same as section 5.5, except the set of (a, b, c, d)
is to express |E*| in the sigmoid function. Since the VE solution is calculated in the time domain,
the converting from |E*| to E(t) is required, assuming the phase angel function φ(ω) may be
unknown.

5.6.1 Forward calculation of E(t) from |E*| and φ
If the mastercurve |E*(ω)| and φ(ω) are obtained from the experiment data, the coefficients
Ei (spring stiffness) of the Generalized Maxwell model can be directly calculated by the inverse
of the coefficient matrix A, as:

E = A−1 * b = ( A'*A) −1 * ( A'*b) ,
where:

[ ]

A = aij

(48)

(ω T )
=
1 + (ω T )
2

and

aij

i

j

2

i

j

[ ] and b = E' (ω ) − E' (∞) =| E * (ω )| *sin(ϕ(ω ))− | E * (∞) | *sin(ϕ(∞))

b = bj

j

j

j

j

Then the relaxation modulus E(t) can be expressed by Ei, as Eq. (4).
However, in the back calculation problem, this procedure is invalid, because it assumes that
in laboratory testing the phase angle φ(ω) is known, but φ(ω) is not known in the case of
backcalculation.

107

One potential solution to this problem is to increase the total number of variables in the
backcalculation procedure from 4 parameters (a, b, c, d) to 8 parameters (a, b, c, d, e, f, g, h).
However, this makes the backcalculation problem more complex and less likely to be solved.
Two other possible solutions are discussed in below sections.

5.6.2 Calculation of E(t) and φ from |E*| by Iteration
The dynamic modulus mastercurve |E*(ω)| can be mathematically expressed by coefficients
of the spring stiffness Ei by Generalized Maxwell model, as shown in Eq. (49).

[E ' (ω )]2 + [E" (ω )]2 = [E * (ω )]2

(49)

(ωTi )2
E ' (ω ) = ∑ Ei mi (ω ) = E0 + ∑ Ei
2
1 + (ωTi )
i =0
i =1
N

where:

,

N

N

N

i =0

i =1

E" (ω ) = ∑ Ei ni (ω ) = ∑ Ei

ωTi
2
1 + (ωTi )

Eq. (49) is a continuous function of frequency, ω, and there are N unknown variables, so the
minimum Root Mean Square Error (RMSE) method can be applied to solve for the N unknown
variables, by discrete frequency ω for M equations. The mathematical expression is:
2
2
2
⎧M ⎡ N
N
⎤ ⎫
⎛
⎞ ⎛
⎞
2
⎪
⎪
Min[F (E )] = Min ⎨∑ ⎢⎜ ∑ mkj Ek ⎟ + ⎜ ∑ nkj Ek ⎟ − E * j ⎥ ⎬
⎠
⎠ ⎝ k =1
⎥ ⎪,
⎪ j =1 ⎢⎝ k =1
⎣
⎦ ⎭
⎩

(ω T )

2

where:

mkj =

j k

1 + (ωTk )

2

and

mk 0 = 1
108

(50)

nkj =

ω jTk
2
1 + (ωTk )

and

nk 0 = 0

By finding the Stationary point of function F(E), a set of M equations can be obtained. In
order to avoid the problem of indeterminate system of equations, one can assume M = N for
simplification. Following the Gauss-Seidel iteration method, or the iteration method of
successive replacement, one can solve for Ek by equation number k in the equation set.
N
N
⎞
⎛
∂F
= 2 K j ⎜ 2mkj ∑ mkj Ek + 2nkj ∑ nkj Ek ⎟ ≡ 0 ,
∂Ek
k =1
k =1
⎝
⎠
2

⎞
⎞ ⎛
⎛
2
K j = ⎜ ∑ mkj Ek ⎟ + ⎜ ∑ nkj Ek ⎟ − E * j
⎠
⎠ ⎝ k =1
⎝ k =1
N

where:

2

(51)

N

The above Eq. (51) translates to a set of cubic polynomial functions of Ei, and
mathematically there are three roots for each equation. It is possible to solve the above set of
equations by iteration, i.e., assuming the seed value of Ei, and then calculating each coefficient
Kj. Thus, the above equation can be simplified to one set of linear equations, and it can be
uniquely identified for new values of Ei. By numerical iteration, a better set of Ei can be obtained
by inversing the coefficient matrix with the minimum RMSE values.
In case of non-invertible matrix, one variable Ei can be solved in the ith equation, assuming
Kj and the rest Ej are known from the seed value. The mathematical equation is expressed in Eq.
(52), and the algorithm is shown below briefly in Fig. 53:

109

N

Ek = −

N

∑

m=1,m≠ k

∑ K (m
j

jk

j =1

N

m jm + n jk n jm )

(

K j m2 + n2
∑
jk
jk
j =1

)

Em

(52)

It is important to note that this procedure is a proposed idea that has not been validated. One
concern is that the cubic polynomial equation is simplified into a set of linear equations, which
may imply that convergence may not be achieved.

110

Input: |E*|
(0 )
Assumed: Ei

i = i +1

Kj by Eq. (51)

Ek(i ) by Eq. (52)
(i)

Update Ei , and calculate

Error (i )

E *(i ) (ω ) , φ (i ) (ω )

(

)

⎡ abs E *i(i ) (w j ) − E *i (w j ) ⎤
= ∑⎢
⎥
E *i (w j )
⎢
⎥
⎣
⎦

Error<ε

2

N

Y
Coefficients for springs (Ei) in generalized Maxwell model

Output: Relaxation Modulus E (t), φ (ω )

Figure 53: Algorithm for Calculation of E(t) and φ from |E*| by iteration

After programming in MATLAB, it was found that the convergence of this method is

111

questionable. As the numerical solution approaches the exact solution, Kj value in Eq. (51)
approaches zero, and the difference between numerical solution and exact solution will be
amplified by Eq. (52) as Kj is in the denominator. Mathematically, the error from simplification
is too large to converge in calculation, as the cubic polynomial equation is simplified into a set of
linear equations.

5.6.3 Calculation of E(t) from |E*(ω)| by Iteration φ(ω)
There is no problem in convergence for the first method in Section 5.6.1, while there is
iteration for the second method in Section 5.6.2 without theoretical confliction. Hence, a better
way is to combine the advantage of the above two methods. Another proposed method is
introduced in this section, and the algorithm is shown below briefly in Fig. 54.

112

Input: |E*|
Assumed: φ

(0 )

(ω )
i = i +1

E'(i ) (ω )

Generalized Maxwell Model: E (i )

φ (i+1) (ω )

E"(i ) (ω )

(

)

⎡ abs E *i(i) (wj ) − E *i (wj ) ⎤
(i )
Error = ∑ ⎢
⎥
E *i (wj )
⎢
⎥
⎣
⎦

Error<ε

2

N

Y
Coefficients for springs (Ei) in generalized Maxwell model

Output: Relaxation Modulus E (t), φ (ω )

Figure 54: Algorithm for Calculation of E(t) from |E*| by iteration φ

After many cases of running, results for the identification of the phase angle φ are
questionable, however, the magnitude of E(t) works well at low reduced time. The dynamic
modulus mastercurve |E*| of one example before and after the conversion is shown in Fig. 55.
113

Magnitude of Dynamic
Modulus (Pa)

1.E+11
1.E+10
Input |E*|
Pred |E*|

1.E+09
1.E+08
1.E+07

1.E+06
1.E-10 1.E-06 1.E-02 1.E+02 1.E+06 1.E+10
Log(Reduced Frequency)
Figure 55: |E*| input and calculated after converting to E(t) for one AC mixture

Based on above discussion, there are only 4 parameters (a, b, c, d) required in
backcalculation for |E*| for the AC layer, and the phase angle shift function φ(ω) with parameters
(e, f, g, h) are intermediately calculated in the program.

5.6.4 Flowchart of VE backcalculation of |E*|
Similar considerations in Section 5.5 are required in this section, but they are ignored here
for brevity. The flow chart of VE backcalculation with AC material property expressed as |E*| is
shown in Fig. 56.

114

Input: Pavement structure profile,
FWD sensor profile, F(t), Wi(t)
Time delay ∆ti, deflection history Wi (tj)

Estimate subgrade modulus via k-value

Seed value:

i=i+1

{X }

AC layer Relaxation Modulus E

(i)

Boundary Check

(t), φ(ω)

Y
Deflection basin W (i)(t ) from VE solution
i
j

⎡
⎛ (i) ⎛ ⎞
⎛ ⎞⎞⎤
⎢ abs ⎜ w i ⎜ t j ⎟ − w i ⎜ t j ⎟ ⎟ ⎥
(i )
⎝ ⎠
⎝ ⎠⎠⎥
⎝
Error = ∑ ⎢
⎢
⎥
w ⎛t ⎞
i ⎜ j⎟
⎢
⎥
⎝ ⎠
⎣
⎦

(i)

Error <ε

N

i<NMax
2

Avoid Local Minimum

Increment by Newton’s Method

Y
New Seed value:

{X }

MATLAB “fminsearch”

Output:

{X} , φ(ω)

Figure 56: Algorithm of backcalculation procedure of |E*|

115

N

5.7 One Numerical Example of Backcalculation Result
A hypothetical pavement profile is assumed, and the response of the pavement under FWD
test is simulated by layered VE solution in Chapter 4. The estimation of elastic modulus of the
subgrade layer is inappropriate, as the effect of mass inertia cannot be simulated via layered VE
solution. Thus, the elastic modulus of the subgrade layer is directly backcalculated. The result of
the backcalculation program is listed by |E*| or E(t) of AC layer separately.

5.7.1 Result of VE backcalculation of E(t)
The actual pavement profile and the backcalculated result of each layer are listed in Table 33.
The input and backcalculated deflection histories are shown in Fig. 57, and the input and
backcalculated E(t) functions of the AC layer are shown in Fig. 58.
Table 33: Basic information of a three-layer pavement structure and backcalculated result

Physical
Layer
AC
Base
Subgrade

Elastic
Modulus
(MPa (ksi))
E*(f)/ E(t)
100 (14.5)
35.0 (5.1)

Poisson’s
Ratio
0.35
0.35
0.45

Thickness
(m (in))
0.1 (4)
0.3 (12)
Infinity

116

E(t)
|E*|
Backcalculation Backcalculation
(MPa (ksi))
(MPa (ksi))
Fig. 56
Fig. 58
98.6 (14.3)
106.7 (15.5)
35.0 (5.1)
34.9 (5.1)

Deflection history simulation (E(t))

deflection (inch)

0.06

sensor 1
sensor 3
sensor 5

0.05
0.04

sensor 2
sensor 4
sensor 6

0.03
0.02
0.01
0.00
0

10

20

30

40

50

time (ms)

Modulus (Pa)

Figure 57: Input and backcalculated time-histories of E(t) in the numerical FWD test

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+10
Backcalculated
Actual
1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

1.E-03

1.E-01

1.E+01

1.E+03

Time (s)
Figure 58: Input and backcalculated E(t) for the numerical FWD test

117

Table 33 shows that the error of the elastic modulus for base layer is around 1.5%, while the
error for subgrade layer is 0.0%. Fig. 57 illustrates that the simulated deflection history matches
very well with the input deflection history, thus, the program converges in calculation. Fig.58
illustrates that the backcalculated E(t) does not match very well with the actual input E(t)
function. However, this result is expected because the critical range in the E(t) function is 0.0001
- 0.1 second, which is used in both the forward and backcalculation program. The result of the
long-time response part of E(t) can be improved if the AC mixture design is given, based on
Witczak’s regression equation. Therefore, this proposed backcalculation algorithm works
reasonably well in this problem.

5.7.2 Result of VE backcalculation of |E*|
The actual pavement profile and the backcalculated result of each layer are listed in Table 33
in above section. The input and backcalculated deflection history is shown in Fig. 59, and the
input and backcalculated E(t) of the AC layer is shown in Fig. 60.

118

Deflection history simulation (|E*|)

deflection (inch)

0.06

sensor 1
sensor 3
sensor 5

0.05

sensor 2
sensor 4
sensor 6

0.04
0.03
0.02
0.01
0.00
0

10

20

30

40

50

time (ms)
Figure 59: Input and backcalculated time-histories of |E*| in the numerical FWD test

|E*| input vs. backcalculated
Backcalculated
Actual

Dynamic Modulus (MPa)

1.E+05
1.E+04
1.E+03
1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Frequency (Hz)
Figure 60: Input and backcalculated |E*| for the numerical FWD test

119

Table 33 shows that the error of the elastic modulus via |E*| backcalculation for base layer is
around 7%, while the error for subgrade layer is around 0.3%. Fig. 59 illustrates that the
simulated deflection history matches very well with the input deflection history, thus, the
program converges in calculation. Fig. 60 illustrates that the backcalculated |E*| does not match
very well with the actual input |E*| function. Similarly, this result is expected because the
correspondingly critical range in the |E*| function is 10 - 10000 Hz. The result of the
low-frequency part of |E*| can be improved by Witczak’s regression equation if the AC mixture
design is given.
Therefore, this proposed backcalculation algorithm works reasonably well in this problem.

5.8 Summary
The output variables of the backcalculation algorithm are identified, and the
corresponding boundaries are found.
The theoretical background is briefly investigated. Newton’s method is chosen to scan
all possible solutions, random function is used to avoid the problem of local minimum,
and the results from MATLAB internal function “fminsearch” are chosen as the final
results.
A step-by-step procedure for backcalculation algorithm is provided in this chapter, and
both cases of |E*| and E(t) backcalculation for the AC mixture are independently
investigated. The accuracy of the proposed algorithm is numerically checked.

120

Chapter 6. Verification and Sensitivity Analysis
6.1 Introduction
So far, only hypothetical pavement responses were used in the backcalculation
algorithm (Chapter 5). In this chapter, more numerical simulations from both VE forward
solution and dynamic solutions, SAPSI and LAMDA, are used as the input for the
backcalculation algorithm, and the accuracy of the program is checked against the actual
modulus values. One example (Texas site) of a three-layer pavement is used to check the
limitations of the new backcalculation algorithm. Additionally, field FWD test deflection
histories are input for the backcalculation program, and the accuracy of the backcalculation
results are checked against the lab measured values. Finally, a sensitivity analysis is conducted
for each input variable, and only the viscoelastic parameters of the AC layer are backcalculated.
The critical variables are found from sensitivity analysis.

6.2 Verification of the VE Backcalculation Algorithm
In this section, the results from backcalculation algorithm of |E*| and E(t) for the AC layer
are proposed together in each example. The parameters in the sigmoid function for the AC layer,
and elastic modulus of base and subgrade layers are the final output results.
The verification part includes three separate methods. First, four VE forward simulation
results are input in the backcalculation program to check the accuracy of the algorithm. Second,
four actual pavements are investigated (Ji, 2005), and the FWD tests are simulated by the

121

dynamic solution in SAPSI and LAMDA. The simulated deflection histories, as well as the layer
thickness and sensor locations, are entered in the backcalculation program. Third, three field
FWD test deflection histories are used to check the potential usage of the algorithm in the future.
The FWD test is usually conducted at a specific temperature, and the backcalculation is done
for the same temperature. There is no need for time temperature superposition, and the actual
time scale is the reference time scale. In the field FWD test, the typically sampling time interval
is 0.4ms, and the total recording time is 60 ms. The available range from FWD test should be
between 0.4 ms to 60 ms, i.e., the valid range of the E(t) for the AC layer is only the proportion
between 0.4 ms to 60 ms. Although the E(t) function is characterized by sigmoidal function,
theoretically only a portion of E(t) can be obtained from backcalculation, shown in Fig. 61.

Relaxiation Modulus E(t)
Relaxation Modulus
(MPa)

1.E+05
E(t)

1.E+04
1.E+03
1.E+02
Available Range

1.E+01

1.E+00
1.E- 1.E- 1.E- 1.E- 1.E- 1.E+ 1.E+ 1.E+ 1.E+ 1.E+ 1.E+
05 04 03 02 01 00 01 02 03 04 05
Time (s)
Figure 61: Available range in E(t) function from field FWD test backcalculation

Approximately, the available range of |E*| in the frequency domain can be identified by

122

f=1/(2*t), i.e., the valid range of the |E*| for the AC layer is only the proportion between 8 Hz to
1250 Hz, if it is backcalculated from field FWD test, shown in Fig. 62.

Dynamic Modulus (MPa)

Dynamic Modulus |E*|
1.E+05

|E*|

1.E+04
1.E+03
1.E+02
1.E+01

Available Range

1.E+00
1.E- 1.E- 1.E- 1.E- 1.E- 1.E+0 1.E+0 1.E+0 1.E+0 1.E+0
05
04
03
02
01
0
1
2
3
4
Time (s)
Figure 62: Useful range in |E*| function from field FWD test backcalculation

The modification of the FWD test equipment may be required if the useful range needs to be
enlarged. In this dissertation, the typical FWD test equipment is utilized in the backcalculation
algorithm, and proportion of E(t) or |E*| in the useful range is concerned.

6.2.1 VE Backcalculation using Time Histories from VE Solution
Four numerical examples of FWD test deflection histories are simulated by the layered VE
solution, which was discussed in Chapter 4. Then, the deflection histories from sensors and
pavement structural information without modulus are entered into the VE backcalculation
program to check the accuracy of the backcalculation results.
123

6.2.1.1 Numerical Example 1
The pavement response under FWD dynamic load is simulated by VE forward solution, and
the deflection time-history of each sensor is used as input to the VE backcalculation program. In
order to compare the accuracy of this VE program against common elastic backcalculation
programs, the deflection basin is used to backcalculate the elastic modulus of each physical layer
by MODCOMP5. The actual pavement profile and the backcalculated result of each layer are
listed in Table 34. The input and VE backcalculated deflection histories using E(t) are shown in
Fig. 63, the input and VE backcalculated deflection histories using |E*| are shown in Fig. 65; the
input and backcalculated E(t) of the AC layer are shown in Fig. 64, and the input and
backcalculated |E*| of the AC layer are shown in Fig. 66.
Table 34 ： Basic information of pavement structure and backcalculated result for
numerical example 1

Elastic
MODCOMP5
E(t)
|E*|
Poisson Thickness
Modulus
Backcalculation Backcalculation Backcalculation
Ratio (m (in))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
AC
E*(f)/ E(t) 0.35 0.152 (6.0) 8067 (1170.0)
Fig. 64
Fig. 66
Base 206.8 (30.0) 0.30 0.670 (24) 215.1 (31.2)
207.7 (30.3)
236.6 (34.3)
Subgrade 34.5 (5.0) 0.45
Infinity
34.3 (5.0)
34.5 (5.0)
33.7 (4.9)
Physical
Layer

Table 34 shows the elastic modulus of base and subgrade layer are almost the same from VE
backcalculation and MODCOMP5, although the VE backcalculation program using E(t) may
give slightly worse results for the base elastic modulus than that using |E*|. However,
MODCOMP5 can only backcalculate one value for the modulus of the AC layer, which can be
considered as an effective value, but cannot be taken as a fundamental property of the AC
material. The VE backcalculation program, on the other hand, can overcome this problem and

124

gives the AC material property as a function of time or frequency.

Deflection history simulation (E(t))

deflection (inch)

2.5E-02

sensor
sensor
sensor
Sensor

2.0E-02

1
3
5
7

sensor 2
sensor 4
sensor 6

1.5E-02
1.0E-02
5.0E-03
0.0E+00
0

10

20

30

40

50

time (ms)
Figure 63: Input and simulated time-histories by backcalculation of E(t) in example 1
(Solid lines represent the input deflection history, symbols represent the VE simulation)
Figure 64：Input and backcalculated E(t) in example 1

Modulus (Pa)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+11
Backcalculated
Actual
1.E+10
1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

1.E-03

1.E-01
Time (s)

125

1.E+01

1.E+03

Figures 63 shows an excellent match between the actual and predicted deflection time
histories using E(t) for every sensor. Figure 64 shows that the backcalculated function E(t) is
very close to the input function only in the concerned time range around 0.0001 – 0.1 s.

Deflection history simulation (|E*|)

deflection (inch)

2.5E-02

sensor
sensor
sensor
Sensor

2.0E-02

1
3
5
7

sensor 2
sensor 4
sensor 6

1.5E-02
1.0E-02
5.0E-03
0.0E+00
0

10

20

30

40

50

time (ms)
Figure 65: Input and simulated time-histories by backcalculation of |E*| in example 1
(Solid lines represent the input deflection history, symbols represent the VE simulation)

126

|E*| input vs. backcalculated
Dynamic Modulus (MPa)

1.E+05
1.E+04
Backcalculated
Actual
E(t) Converted

1.E+03
1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Frequency (Hz)
Figure 66: Input and backcalculated |E*| in example 1

Figure 65 shows a very good match between actual and predicted deflection time histories
using |E*| for every sensor. Figure 66 shows that the backcalculated function of |E*| for the AC
layer is very close to the input function in the entire frequency range. It is not expected that the
predicted function matches well with the actual |E*| function for every frequency, since only a
portion of the |E*| function is used in both forward and backcalculation solution.
As discussed in Chapter 2, the |E*| can be mathematically calculated from E(t) function. If
the E(t) function is backcalculated as Figure 64, the corresponding |E*| curve is plotted in Figure
66. The coverted |E*| matches very well with the actual |E*| function in the concerned frequency
range, i.e., 10 to 10000 Hz.
127

Phase angle shift (o)

Phase Angle shift function: Backcalculated Vs.
Actual
Backcalculated
Actual
E(t) Converted

50
40
30
20
10
0
1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Frequency (Hz)
Figure 67: Input and backcalculated φ in example 1

Figure 67 shows that the backcalculated phase angle φ for the AC material matches fairly to
the input function in the concerned range. Another predicted φ was obtained from the
interconversion of E(t), as discussed in Chapter 2. Although the predicted E(t) matches very well
to the input function in the concerned range (Figure 67), the converted φ does not match the
input function. In general, it seems φ cannot be well predicted in either E(t) backcalculation, or
|E*| backcalculation, although it works in above example. Therefore, the backcalculation of φ is
ignored in the future examples.
In this example, the VE backcalculation program works well for both |E*| and E(t)
backcalculation, and the match between backcalculated and actual deflection time-histories is
very good. The elastic modulus for the subgrade layer is very accurate, with an error of only 1%;

128

the function E(t) and |E*| for AC layer match fairly well with the input functions in the VE
forward calculation; the elastic modulus of the base layer may have an error of up to 10%.

6.2.1.2 Numerical Example 2
Similarly to the first example, the pavement response under FWD dynamic load is simulated
by the VE forward solution, and the deflection time-history of each sensor is used as input to the
VE backcalculation program. However, in this example, a stiff cement treated base (CTB) is
included in the pavement structure. In order to compare the accuracy of the VE program against
common elastic backcalculation program, the deflection basin is used to backcalculate the elastic
modulus of each physical layer by MODCOMP5. The actual pavement profile and the
backcalculated results of each layer are listed in Table 35. The input and VE backcalculated
deflection histories of E(t) are shown in Fig. 68; the input and VE backcalculated deflection
histories of |E*| are shown in Fig. 70; the input and backcalculated E(t) of the AC layer are
shown in Fig. 69; and the input and backcalculated |E*| of the AC layer are shown in Fig. 71.
Table 35 ： Basic information of pavement structure and backcalculated result for
numerical example 2

Elastic
MODCOMP5
E(t)
|E*|
Poisson Thickness
Modulus
Backcalculation Backcalculation Backcalculation
Ratio (m (in))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
AC
Fig. 69
Fig. 71
E*(f)/ E(t) 0.35 0.102 (4.0) 10618 (1540)
Base
5000 (725) 0.20 0.152 (6.0) 4165 (604)
5104 (740)
5153 (747)
Subgrade 69.0 (10.0) 0.45
Infinity
69.6 (10.1)
72.1 (10.5)
70.8 (10.3)
Physical
Layer

Table 35 shows the VE backcalculated modulus of base and subgrade layers are almost the
same as the actual value, and MODCOMP5 gives slightly worse result for the CTB layer, with an
error of about 20%. Also, MODCOMP5 can only backcalculate one value for the modulus of the
129

AC layer, which can be considered as an effective value, but cannot be taken as a fundamental
property of the AC material. The VE backcalculation program, on the other hand, can overcome
this problem and gives the AC material property as a function of time or frequency.

Deflection history simulation (E(t))

deflection (inch)

1.6E-02

sensor
sensor
sensor
Sensor

1.4E-02
1.2E-02

1
3
5
7

sensor 2
sensor 4
sensor 6

1.0E-02
8.0E-03
6.0E-03
4.0E-03
2.0E-03
0.0E+00
0

10

20

30

40

50

time (ms)
Figure 68: Input and simulated time-histories by backcalculation of E(t) in example 2
(Solid lines represent the input deflection history, symbols represent the VE simulation)

130

Modulus (Pa)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+11
Backcalculated
Actual
1.E+10
1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

1.E-03

1.E-01

1.E+01

1.E+03

Time (s)
Figure 69：Input and backcalculated of E(t) in example 2

Fig. 68 shows a very good match between actual and predicted deflection time histories
using E(t) for every sensor. Fig. 69 shows that the backcalculated function E(t) is very close to
the input function in the concerned time range, approximately 0.0001 – 0.1 s.

131

Deflection history simulation (|E*|)

deflection (inch)

1.6E-02

sensor
sensor
sensor
Sensor

1.4E-02
1.2E-02

1
3
5
7

sensor 2
sensor 4
sensor 6

1.0E-02
8.0E-03
6.0E-03
4.0E-03
2.0E-03
0.0E+00
0

10

20

30

40

50

time (ms)
Figure 70: Input and simulated time-histories by backcalculation of |E*| in example 2
(Solid lines represent the input deflection history, symbols represent the VE simulation)

|E*| input vs. backcalculated

Dynamic Modulus (MPa)

1.E+05
1.E+04
1.E+03
Backcalculated
Actual

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02

1.E+00

1.E+02

Frequency (Hz)
Figure 71: Input and backcalculated of |E*| in example 2

132

1.E+04

Fig. 70 shows a very good match between actual and predicted deflection time histories
using |E*| for every sensor. Fig. 71 illustrates that the backcalculated function |E*| for the AC
layer is very close to the input function in the concerned frequency range, 10 – 10000 Hz. Here
the lower frequency part of |E*| is quite different between input and backcalculated results.
In summary, the VE backcalculation program works well for both |E*| and E(t)
backcalculation, and the match between actual and backcalculated deflection time-histories is
very good. The elastic modulus for the CTB base and subgrade layer is very accurate, with an
error of 3%, and the functions E(t) and |E*| for the AC layer match fairly well with the input
function within certain time and frequency ranges.

6.2.1.3 Numerical Example 3
In this example, a weak pavement section is used. The actual pavement profile and the
backcalculated results of each layer are listed in Table 36. The input and VE backcalculated
deflection histories of E(t) are shown in Fig. 72; the input and VE backcalculated deflection
histories of |E*| are shown in Fig. 74; the input and backcalculated E(t) of the AC layer are
shown in Fig. 73, and the input and backcalculated |E*| of the AC layer are shown in Fig. 75.
Table 36 ： Basic information of pavement structure and backcalculated result for
numerical example 3

Elastic
MODCOMP5
E(t)
|E*|
Poisson Thickness
Modulus
Backcalculation Backcalculation Backcalculation
Ratio (m (in))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
AC
Fig. 73
Fig. 75
E*(f)/ E(t) 0.35 0.102 (4.0) 27028 (3920)
Base
10.8 (1.57) 0.35 0.204 (8.0)
124 (18.0)
10.45 (1.52)
10.9 (1.58)
Subgrade 10.0 (1.45) 0.35
Infinity
13.2 (1.92)
10.1 (1.47)
9.76 (1.42)
Physical
Layer

133

Table 36 shows the VE backcalculated elastic modulus of base and subgrade layer are almost
the same as the actual values, but MODCOMP5 gives quite different result for the base layer.
Moreover, only one modulus value is given in MODCOMP5 for the AC layer. The VE
backcalculation program can overcome this problem and gives the AC modulus as a function of
time or frequency.

Deflection history simulation (E(t))

deflection (inch)

7.E-02

sensor 1
sensor 3
sensor 5

6.E-02
5.E-02

sensor 2
sensor 4
sensor 6

4.E-02
3.E-02
2.E-02
1.E-02
0.E+00
0

10

20

30

40

50

time (ms)
Figure 72: Input and simulated time-histories by backcalculation of E(t) in example 3
(Solid lines represent the input deflection history, symbols represent the VE simulation)

Fig. 72 shows a very good match between actual and predicted deflection time histories
using E(t) for every sensor, with the exception of sensor 1 in the rest period. This is why the
backcalculated function E(t) is very close to the input function only in the concerned time range
0.0001 – 0.1 s (Fig. 73).

134

Modulus (Pa)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+11
Backcalculated
Actual
1.E+10
1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

1.E-03

1.E-01

1.E+01

1.E+03

Time (s)
Figure 73：Input and backcalculated E(t) in example 3

Deflection history simulation (|E*|)

deflection (inch)

7.E-02

sensor 1
sensor 3
sensor 5

6.E-02
5.E-02

sensor 2
sensor 4
sensor 6

4.E-02
3.E-02
2.E-02
1.E-02
0.E+00
0

10

20

30

40

50

time (ms)
Figure 74: Input and simulated time-histories by backcalculation of |E*| in example 3
(Solid lines represent the input deflection history, symbols represent the VE simulation)

135

|E*| input vs. backcalculated

Dynamic Modulus (MPa)

1.E+05
1.E+04
1.E+03

Backcalculated
Actual

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Frequency (Hz)
Figure 75: Input and backcalculated |E*| in example 3

Fig. 74 and 75 show good agreement between deflection time histories and |E*| curves. The
simulated time-history for sensor 1 is almost the same as the input deflection history, so the
backcalculated |E*| function for the AC layer is very close to the input function in the concerned
frequency range 10 – 10000 Hz.
In summary, the VE backcalculation program works well for both |E*| and E(t)
backcalculation, and the match between actual and backcalculated deflection time-histories is
very good. The backcalculated elastic modulus for the base and subgrade layer is very accurate,
with an error of 3%, and the function E(t) and |E*| for AC layer match fairly well with the actual
function for limited ranges in time and frequency. Since the moduli of the base and subgrade
layer are very similar, MODCOMP5 does not work well in this case.
136

6.2.1.4 Numerical Example 4
In this fourth example, the same pavement section as the one in the previous example is used,
except that a different AC mixture is used. The actual pavement profile and the backcalculated
result of each layer are listed in Table 37. The input and VE backcalculated deflection histories of
E(t) are shown in Fig. 76; the input and VE backcalculated deflection histories of |E*| are shown
in Fig. 78; the input and backcalculated E(t) of the AC layer are shown in Fig. 77; and the input
and backcalculated |E*| of the AC layer are shown in Fig. 79.
Table 37 ： Basic information of pavement structure and backcalculated result for
numerical example 4

Elastic
MODCOMP5
E(t)
|E*|
Poisson Thickness
Modulus
Backcalculation Backcalculation Backcalculation
Ratio (m (in))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
AC
Fig. 77
Fig. 79
E*(f)/ E(t) 0.35 0.102 (4.0) 9515 (1380)
Base
10.8 (1.57) 0.35 0.204 (8.0) 103.4 (15.0)
10.5 (1.52)
10.6 (1.54)
Subgrade 10.0 (1.45) 0.35
Infinity
9.45 (1.37)
9.9 (1.44)
9.80 (1.42)
Physical
Layer

Table 37 shows the VE backcalculated elastic modulus of base and subgrade layer are almost
the same from the actual value, but MODCOMP5 gives quite different result for the base layer,
due to the fact that the base modulus is almost the same as the subgrade layer. Moreover, only
one modulus value is given in MODCOMP5 for AC layer. The VE backcalculation program can
overcome this problem and gives the AC modulus as a function of time or frequency.

137

Deflection history simulation (E(t))

deflection (inch)

1.2E-01

sensor 1
sensor 3
sensor 5

1.0E-01
8.0E-02

sensor 2
sensor 4
sensor 6

6.0E-02
4.0E-02
2.0E-02
0.0E+00
0

10

20
30
time (ms)

40

50

Figure 76: Input and simulated time-histories by backcalculation of E(t) in example 4

Modulus (Pa)

(Solid lines represent the input deflection history, symbols represent the VE simulation)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+10
Backcalculated
Actual
1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

Figure 77：Input and backcalculated E(t) in example 4

138

1.E+03

Fig. 76 and 77 show comparisons of backcalculated and actual deflection time histories and
E(t) curve. The backcalculated E(t) function is very close to the input function only in the
concerned time range, 0.0001 – 0.01 s.

Deflection history simulation (|E*|)
deflection (inch)

1.2E-01

sensor 1
sensor 3
sensor 5

1.0E-01

sensor 2
sensor 4
sensor 6

8.0E-02
6.0E-02
4.0E-02
2.0E-02
0.0E+00
0

10

20
30
time (ms)

40

50

Figure 78: Input and simulated time-histories by backcalculation of |E*| in example 4
(Solid lines represent the input deflection history, symbols represent the VE simulation)

139

|E*| input vs. backcalculated
Dynamic Modulus (MPa)

1.E+04
1.E+03
Backcalculated
Actual

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 79: Input and backcalculated |E*| in example 4

Fig. 78 and 79 illustrate good agreements between backcalculated and actual deflection time
histories and |E*| curves. The simulated time-history for sensor 1 is close to the input deflection
history, so the backcalculated function for the AC layer is reasonably close to the input function
in the concerned frequency range 10 – 10000 Hz, although larger errors exist in other frequency
ranges.
In summary, another AC mixture is used in this numerical example, and the backcalculation
results are still reasonable. The elastic moduli for the base and subgrade layer are very accurate,
with an error of 3%, and the functions E(t) and |E*| for AC layer match fairly well with the input
functions in the concerned time and frequency ranges. Since the modulus of the base and
subgrade layer is very similar, MODCOMP5 does not work well in this case.

140

6.2.1.5 Discussion
The above four numerical examples indicate that the backcalculation algorithm works well,
even when there is a stiff base. Since the forward simulation is from VE solution, the dynamic
effect of FWD test cannot be shown, and the elastic modulus of the subgrade layer is obtained
directly from backcalculation, not using the ratio of shear and Rayleigh wave velocity discussed
in Chapter 3.

6.2.2 VE Backcalculation using Time Histories from Dynamic
Solutions
In the previous section, the layered VE solution was used for both forward and
backcalculation algorithm. This allowed for the verification of the inverse solution. In this
section, we use completely different solution to calculate the deflection time histories under
FWD test. Four numerical examples of FWD test deflection histories are simulated by dynamic
solution SAPSI or LAMDA. Then, the deflection histories from sensors and pavement structural
information without modulus are entered into VE backcalculation program to check the accuracy
of the backcalculation results. One of the examples includes a site of shallow stiff layer
underneath pavement.

6.2.2.1 Example 5: Michigan Site
The pavement response under FWD dynamic load is simulated by SAPSI, and the deflection
time-history of each sensor is the input to the VE backcalculation program. In order to compare
the accuracy of this VE program against common elastic backcalculation program, the deflection
basin is used to backcalculate the elastic modulus of each physical layer by MODCOMP5. The
141

actual pavement profile and the backcalculated result of each layer are listed in Table 38. The
input and VE backcalculated deflection histories using E(t) are shown in Fig. 80; the input and
VE backcalculated deflection histories using |E*| are shown in Fig. 82; the input and
backcalculated E(t) of the AC layer are shown in Fig. 81; and the input and backcalculated |E*|
of the AC layer are shown in Fig. 83.
Table 38：Basic information of Michigan pavement structure and backcalculated results

Elastic
MODCOMP5
E(t)
|E*|
Poisson Thickness
Modulus
Backcalculation Backcalculation Backcalculation
Ratio (m (in))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
AC
Fig. 81
Fig. 83
E*(f)/ E(t) 0.35 0.178 (7.0) 4240 (614.9)
Base 206.8 (30.0) 0.35 1.321 (52) 206.2 (29.9)
208.6 (30.3)
208.7 (30.3)
Subgrade 310.3 (45.0) 0.45
Infinity
286.1 (41.5)
313.0 (45.4)
312.4 (45.3)
Physical
Layer

Table 38 shows the elastic modulus of base and subgrade layers are almost the same from
VE backcalculation and MODCOMP5, although the VE backcalculation program gives slightly
better results. However, only one modulus value is given in MODCOMP5 for the AC layer. The
VE backcalculation program can overcome this problem and gives the AC modulus as a function
of time or frequency.

142

Deflection history simulation (E(t))

deflection (inch)

1.2E-02

sensor
sensor
sensor
Sensor

1.0E-02
8.0E-03

1
3
5
7

sensor
sensor
sensor
Sensor

2
4
6
8

6.0E-03
4.0E-03
2.0E-03
0.0E+00
0

10

20
30
time (ms)

40

50

Figure 80: Input and simulated time-histories by backcalculation of E(t) in Michigan site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

Relaxation Modulus Comparison: Backcalculated
1.E+11
Vs. Input
Backcalculated
Actual

Modulus (Pa)

1.E+10
1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

Figure 81：Input and backcalculated E(t) in Michigan site

143

1.E+03

Fig. 80 and 81 show good agreement between backcalculated and actual deflection time
histories and E(t) curves. The difference between input and simulated time-history for sensor 1 is
not significant, so the backcalculated and input E(t) mach well (Fig. 81), especially in the

deflection (inch)

relevant time range 0.0001 – 0.1 s.

Deflection history simulation (|E*|)
1.2E-02
sensor 1
sensor 3
1.0E-02
sensor 5
Sensor 7
8.0E-03

sensor
sensor
sensor
Sensor

2
4
6
8

6.0E-03
4.0E-03
2.0E-03
0.0E+00
-2.0E-03
0

10

20
30
time (ms)

40

50

Figure 82: Input and simulated time-histories by backcalculation of |E*| in Michigan site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

144

|E*| input vs. backcalculated
Dynamic Modulus (MPa)

1.E+05
1.E+04
1.E+03
Backcalculated
Actual

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 83: Input and backcalculated |E*| in Michigan site

Fig. 82 and 83 illustrate good agreement between backcalculated and actual deflection time
histories and |E*| curves. The difference between input and simulated time-history for sensor 1 is
good enough, so the backcalculated |E*| function for the AC layer is very close to the input
function in the concerned frequency range 10 to 10000 Hz.
In summary, the VE backcalculation program works well for both |E*| and E(t)
backcalculation, and the backcalculated results are to input deflection time-histories. The elastic
moduli for the base and subgrade layer are very accurate, with an error of 1%, while the E(t) and
|E*| functions for AC layer match fairly well with the input function in SAPSI.

6.2.2.2 Example 6: Texas Site
A similar procedure is done for one pavement profile in Texas site, where there is a shallow
stiff layer underneath. The actual pavement profile and the backcalculated result of each layer are
145

listed in Table 39. The input and VE backcalculated deflection histories using E(t) are shown in
Fig. 84; the input and VE backcalculated deflection histories using |E*| are shown in Fig. 86; the
input and backcalculated E(t) of the AC layer are shown in Fig. 85; and the input and
backcalculated |E*| of the AC layer are shown in Fig. 87.
Table 39：Basic information of Texas pavement structure and backcalculated results

Elastic
MODCOMP5
E(t)
|E*|
Poisson Thickness
Modulus
Backcalculation Backcalculation Backcalculation
Ratio (m (in))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
AC
Fig. 85
Fig. 87
E*(f)/ E(t) 0.35 0.203 (8.0) 4206 (610.0)
Base 172.4 (25.0) 0.45 1.702 (67.0) 152.4 (22.1)
184.6 (27.0)
197.3 (28.6)
Subgrade 861.8 (125) 0.25
Infinity
500.6 (72.6)
903.9 (131.1)
887.7 (128.8)
Physical
Layer

Table 39 shows the elastic modulus of base and subgrade layers are almost the same from
VE backcalculation and MODCOMP5, although the VE backcalculation program gives slightly
better results. However, only one modulus value is given in MODCOMP5 for the AC layer. The
VE backcalculation program can overcome this problem and gives the AC modulus as a function
of time or frequency.

146

Deflection history simulation (E(t))

deflection (inch)

1.2E-02

sensor
sensor
sensor
Sensor

1.0E-02
8.0E-03

1
3
5
7

sensor
sensor
sensor
Sensor

2
4
6
8

6.0E-03
4.0E-03
2.0E-03
0.0E+00
0

10

20
30
time (ms)

40

50

Figure 84：Input and simulated time-histories by backcalculation of E(t) in Texas site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

Modulus (Pa)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+10
Backcalculated
Actual

1.E+09

1.E+08
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

Figure 85：Input and backcalculated E(t) in Texas site

147

1.E+03

Fig. 84 and 85 illustrate that the backcalculated and actual deflection time histories and E(t)
curves do not match well. The reason of this discrepancy is that the stiff layer traps the wave
energy within the pavement structure, and thus it cannot be simulated by the layered VE solution.

Deflection history input and simulation (|E*|)

deflection (inch)

1.2E-02

sensor
sensor
sensor
Sensor

1.0E-02
8.0E-03

1
3
5
7

sensor
sensor
sensor
Sensor

2
4
6
8

6.0E-03
4.0E-03
2.0E-03
0.0E+00
0

10

20
30
time (ms)

40

50

Figure 86：Input and simulated time-histories by backcalculation of |E*| in Texas site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

148

|E*| input vs. backcalculated
Dynamic Modulus (MPa)

1.E+04
1.E+03
Backcalculated

1.E+02

Actual

1.E+01
1.E+00
1.E-06

1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Frequency (Hz)
Figure 87：Input and backcalculated |E*| in Texas site

Fig. 86 and 87 show that the backcalculated and actual deflection time histories and |E*|
curves do not match well, especially for the farther sensors. The difference between input and
simulated time-history for sensor 1 is not significant in the loading region, but the differences
between input and simulated time-histories for other sensors are very different. Even the base
and subgrade elastic moduli are very similar, the backcalculation program is inappropriate in this
example.
In summary, the VE backcalculation program does not work well, because there is shallow
shift layer underneath. As discussed in Chapter 4, the VE forward solution cannot simulate the
pavement response under FWD test if shallow bedrock exists. It is expected that the
backcalculation program does not converge to the right value. Despites these issues, Table 39
149

shows the backcalculation algorithm converges to the right solution for the base and subgrade
layers, with an error of about 3%.

6.2.2.3 Example 7: Florence Site
A similar procedure is done for one pavement structure in Florence. The actual pavement
profile and the backcalculated result of each layer are listed in Table 40. The input and VE
backcalculated deflection time-histories using E(t) are shown in Fig. 88; the input and VE
backcalculated deflection time-histories using |E*| are shown in Fig. 90; the input and
backcalculated E(t) of the AC layer are shown in Fig. 89; and the input and backcalculated |E*|
of the AC layer are shown in Fig. 91.
Table 40: Basic information of a Florence pavement structure and backcalculated result

Elastic
MODCOMP5
E(t)
|E*|
Poisson Thickness
Modulus
Backcalculation Backcalculation Backcalculation
Ratio (m (in))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
AC
Fig. 89
Fig. 91
E*(f)/ E(t) 0.35 0.102 (4.0) 7860.0 (1140)
Base 275.8 (40.0) 0.35 2.438 (96.0) 271.0 (39.3)
280.5 (40.7)
268.9 (39.0)
Subgrade 551.6 (80.0) 0.15
Infinity
385.4 (55.9)
477.2 (69.2)
657.6 (95.4)
Physical
Layer

Table 40 shows the elastic modulus of base layers are almost the same from both VE
backcalculation and MODCOMP5. However, VE backcalculation gives better result for the
modulus of the stabilized subgrade layer. Only one modulus value is given in MODCOMP5 for
the AC layer. The VE backcalculation program can overcome this problem and gives AC
modulus as a function of time or frequency.

150

Deflection history simulation (E(t))

deflection (inch)

1.4E-02

sensor
sensor
sensor
Sensor

1.2E-02
1.0E-02

1
3
5
7

sensor
sensor
sensor
Sensor

2
4
6
8

8.0E-03
6.0E-03
4.0E-03
2.0E-03
0.0E+00
0

10

20
30
time (ms)

40

50

Figure 88: Input and backcalculation simulated time-histories using E(t) in Florence site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

Modulus (Pa)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+11
Backcalculated
Actual

1.E+10

1.E+09

1.E+08
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

Figure 89: Input and backcalculated E(t) in Florence site

151

1.E+03

Fig. 88 and 89 illustrate good agreement between actual and backcalculated deflection time
histories and E(t) curves. The backcalculated and input E(t) do not mach well for long-time part
(Fig. 89). However, the backcalculated function is very close to the input function in the
concerned time range 0.0001 – 0.1 s.

Deflection history simulation (|E*|)

deflection (inch)

1.4E-02

sensor
sensor
sensor
Sensor

1.2E-02
1.0E-02

1
3
5
7

sensor
sensor
sensor
Sensor

2
4
6
8

8.0E-03
6.0E-03
4.0E-03
2.0E-03
0.0E+00
0

10

20
30
time (ms)

40

50

Figure 90: Input and backcalculation simulated time-histories using |E*| in Florence site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

152

|E*| input vs. backcalculated
Dynamic Modulus (MPa)

1.E+05
1.E+04
1.E+03

Backcalculated
Actual

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Frequency (Hz)
Figure 91: Input and backcalculated |E*| in Florence site

Fig. 90 and 91 show reasonable agreement between the backcalculated and actual deflection
time histories and |E*| curves. The difference between input and simulated time-history for
sensor 1 is not significant, and the backcalculated function for the AC layer is very close to the
input function in concerned frequency range 10 – 10000 Hz.
In summary, the VE backcalculation program works well for both |E*| and E(t)
backcalculation, and there is good agreement between input and backcalculated deflection
time-histories. However, the elastic modulus of the subgrade layer is not as accurate, with an
error around 20%, although it is much better than MODCOMP backcalculation result.

153

6.2.2.4 Example 8: Kansas Site
The actual pavement profile and the backcalculated result of each layer are listed in Table 41.
The input and VE backcalculated deflection time-histories using E(t) are shown in Fig. 92; the
input and VE backcalculated deflection time-histories using |E*| are shown in Fig. 94; the input
and backcalculated E(t) of the AC layer are shown in Fig. 93; and the input and backcalculated
|E*| of the AC layer are shown in Fig. 95.
Table 41: Basic information of a Kansas pavement structure and backcalculated result

Elastic
MODCOMP5
E(t)
|E*|
Poisson Thickness
Modulus
Backcalculation Backcalculation Backcalculation
Ratio
(m (in))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
(MPa (ksi))
AC
Fig. 93
Fig. 95
E*(f)/ E(t) 0.30 0.102 (4.0) 5302.1 (769.0)
Base 413.7 (60.0) 0.35 0.356 (14.0) 447.5 (64.9)
469.1 (68.0)
513.6 (74.5)
Subgrade 482.6 (70.0) 0.45
Infinity
524.7 (76.1)
509.3 (73.9)
507.1 (73.5)
Physical
Layer

Table 41 shows the elastic modulus of base and subgrade layers are similar from both VE
backcalculation and MODCOMP5, although the VE backcalculation program gives slightly
better results for the subgrade layer, while MODCOMP5 gives better result for the base layer.
However, only one modulus value is given in MODCOMP5 for the AC layer. The VE
backcalculation program can overcome this problem and gives AC modulus as a function of time
or frequency.

154

Deflection history simulation (E(t))

deflection (inch)

1.E-02

sensor
sensor
sensor
Sensor

8.E-03
6.E-03

1
3
5
7

sensor
sensor
sensor
Sensor

2
4
6
8

4.E-03
2.E-03
0.E+00
0

10

20
30
time (ms)

40

50

Figure 92: Input and simulated time-histories by backcalculation of E(t) in Kansas site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

Modulus (Pa)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+11
Backcalculated
Actual

1.E+10
1.E+09
1.E+08
1.E+07
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

Figure 93: Input and backcalculated E(t) in Kansas site

155

1.E+03

Fig. 92 and 93 show good agreement between actual and backcalculated deflection time
histories and E(t) curves. The backcalculated and input E(t) do not mach well for long-time part
(Fig. 93). However, the backcalculated function is very close to actual input in the concerned
time range 0.0001 – 0.1 s.

Deflection history simulation (|E*|)

deflection (inch)

1.E-02

sensor
sensor
sensor
Sensor

8.E-03

1
3
5
7

sensor
sensor
sensor
Sensor

2
4
6
8

6.E-03
4.E-03
2.E-03
0.E+00
0

10

20
30
time (ms)

40

50

Figure 94: Input and simulated time-histories by backcalculation of |E*| in Kansas site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

156

|E*| input vs. backcalculated
Dynamic Modulus (MPa)

1.E+05
1.E+04
1.E+03
Backcalculated
Actual

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Frequency (Hz)
Figure 95: Input and backcalculated |E*| in Kansas site

Fig. 94 and 95 illustrate good agreement between actual and backcalculated deflection time
histories and |E*| curves. There is almost no difference between input and simulated
time-histories for sensor 1, except for the rest period, and the backcalculated |E*| function for the
AC layer is very close to the input function in the concerned frequency range 10 – 10000 Hz.
In summary, the VE backcalculation program works reasonably well for both |E*| and E(t)
backcalculation, and there is good match between input and backcalculated deflection time
histories. However, the elastic modulus of the base layer using both |E*| and E(t) backcalculation
is not as accurate, with an error around 10 - 20%, although the modulus of the subgrade layer
from VE backcalculation is slightly better than MODCOMP backcalculation result.

157

6.2.2.5 Discussion
In this section, the accuracy of the backcalculation algorithm is investigated by 4 actual
pavement structures with hypothetical FWD test. The backcalculation algorithm does not work
well if there is shallow stiff layer underneath, as shown in the example of Texas site, although the
AC property at the relevant time/frequency range is acceptable.
Since no experimental data available, the base and subgrade moduli are obtained from
static/dynamic backcalculation (Ji, 2005), which may be questionable in cases. For example, the
modulus of the base layer is lower than the modulus of the subgrade layer. This phenomena is
unrealistic. The VE backcalculation algorithm is forced to simulate the static/dynamic
backcalculation result, so the moduli for base and subgrade layers are not accurate in cases, with
error as high as 20%. If the realistic pavement structure is assumed as in the previous section, the
error of the elastic modulus for the base and subgrade should be much smaller.
The |E*| and E(t) backcalculation results are very similar to each other. However, the E(t)
backcalculation is much faster, as there is no need to convert from |E*| to E(t). Therefore, E(t)
backcalculation method is chosen as the VE backcalculation algorithm in future applications.

6.2.3 VE Backcalculation using Field FWD Test Data
Three examples of FWD test deflection histories were obtained from the FHWA’s APT
(Accelerated Pavement Testing) experiments at the Turner Fairbank Highway Research Center
(TFHRC). Deflection histories from sensors and pavement structural information without
modulus were entered into the VE backcalculation program to check the accuracy of the
backcalculation results. Since tests locations are close to each other, the elastic modulus of the
base and subgrade layers are expected to be very similar.
158

6.2.3.1 FWD Test Conducted at L9S4
One record of the field FWD test was found for Lane 9 (AC mixture: SBS 64-40) in the
FHWA’s APT experiments. The pavement structure includes 6 in AC surface, 20 in Crushed
Aggregate Base (CAB) and deep AASHTO A-4 subgrade soil (bedrock is at 25 ft from surface).
The FWD test was conducted at around 10 AM on October 21, 2010, with air temperature of
o

14.5 C. The dynamic modulus |E*| was measured in the lab with the same AC mixture design,
but compacted with standard gyratory compactor at 4% air void, and the temperature shift factor
was also obtained from time-temperature superposition. The lab E(t) was mathmatically obtained
from lab measured |E*|, as discussed in Section 2.4. The actual modulus for each physical layer
is unknown, the pavement structure and the backcalculated result of each layer are listed in Table
42. The input and VE backcalculated deflection time-histories of E(t) are shown in Fig. 96; the
input and backcalculated E(t) of the AC layer are shown in Fig. 97.
Table 42: Basic information of L9S4 pavement structure and backcalculated result

Physical Poisson’s Thickness
MODCOMP5
E(t) Backcalculation
Layer Ratio
(m (in)) Backcalculation (MPa (ksi))
(MPa (ksi))
AC
0.35
0.152 (6.0)
4509 (654)
Fig. 97
Base
0.30
0.508 (20.0)
66.7 (9.7)
69.6 (10.1)
Subgrade
0.45
Infinity
132.4 (19.2)
135.8 (19.7)
Table 42 shows the elastic modulus of base and subgrade layers are similar from both VE
backcalculation and MODCOMP5. However, only one modulus value is given in MODCOMP5
for the AC layer. The VE backcalculation program can overcome this problem and gives AC
modulus as a function of time.

159

deflection (mls)

Deflection history simulation (E(t))
40
35
30
25
20
15
10
5
0
-5

sensor 1
sensor 4
Sensor 7

0

10

20
30
time (ms)

sensor 2
sensor 5
Sensor 8

sensor 3
sensor 6
Sensor 9

40

50

Figure 96: Input and simulated time-histories by backcalculation of E(t) in L9S4 site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

Modulus (Pa)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+10
Backcalculated
Actual

1.E+09

1.E+08

1.E+07
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

Figure 97: Input and backcalculated E(t) in L9S4 site

160

1.E+03

Fig. 96 shows that the predicted deflection time histories match very well with the actual
values for all the sensors. Moreover, the backcalculated and input E(t) does mach well (Fig. 97),
especially during the concerned time range 0.0001 – 0.1 s.

6.2.3.2 FWD Test Conducted at L10S4
A second record of the field FWD test was found for Lane 10 (AC mixture: Air Blown) in
the FHWA’s APT experiments. The pavement structure is the same as for Lane 9 except for the
AC materials. The FWD test was conducted at around 9 AM on October 21, 2010, with air
o

temperature of 12.3 C. The actual modulus for each physical layer is unknown, the pavement
structure and the backcalculated results of each layer are listed in Table 43. The input and VE
backcalculated deflection time-histories of E(t) are shown in Fig. 98; the input and
backcalculated E(t) of the AC layer are shown in Fig. 99.
Table 43: Basic information of L10S4 pavement structure and backcalculated result

Physical Poisson’s Thickness
Layer Ratio
(m (in))
AC
0.35
0.152 (6.0)
Base
0.30
0.508 (20.0)
Subgrade 0.45
Infinity

MODCOMP5
Backcalculation (MPa (ksi))
10687 (1550)
66.9 (9.7)
178.6 (25.9)

E(t) Backcalculation
(MPa (ksi))
Fig. 99
62.4 (9.0)
176.6 (25.6)

Table 43 shows that the elastic moduli of the base and subgrade layers are similar between
VE backcalculation and MODCOMP5. Only one modulus value is given in MODCOMP5 for the
AC layer. The VE backcalculation program can overcome this problem and gives AC modulus as
a function of time.

161

Deflection history simulation (E(t))

deflection (mls)

25

sensor 1
sensor 4
Sensor 7

20

sensor 2
sensor 5
sensor 8

sensor 3
sensor 6
Sensor 9

15
10
5
0
-5
0

5

10

15

20
25
time (ms)

30

35

40

45

Figure 98: Input and simulated time-histories by backcalculation of E(t) in L10S4 site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

Modulus (Pa)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+11
Backcalculated
Actual

1.E+10

1.E+09

1.E+08
1.E-03

1.E-02

1.E-01

1.E+00 1.E+01
Time (s)

1.E+02

Figure 99: Input and backcalculated E(t) in L10S4 site

162

1.E+03

Fig. 98 shows an excellent match between actual and predicted deflection time histories
using E(t) for every sensor. Figure 95 illustrates that the backcalculated E(t) function is very
close to the lab-measured function, especially in the concerned time range 0.0001 – 0.1 s. It
indicates there is a tiny shift between lab measured and backcalculated E(t), mainly because the
temperature information may not be accurate.

6.2.3.3 FWD Test Conducted at L11S4
A third record of the field FWD test was found for Lane 11 (AC mixture: SBS-LG) in the
FHWA’s APT experiments. The pavement structure is the same as that in the previous examples
except for the AC mixture. The FWD test was conducted at around 2 PM on October 22, 2010,
with AC surface temperature of 30.4oC; however, the record air temperature is 14.8oC. The actual
modulus for each physical layer is unknown, the pavement structure and the backcalculated
results of each layer are listed in Table 44. The input and VE backcalculated deflection
time-histories using E(t) are shown in Fig. 100; the input and backcalculated E(t) of the AC layer
are shown in Fig. 101.
Table 44: Basic information of L11S4 pavement structure and backcalculated result

Physical Poisson’s Thickness
MODCOMP5
Layer Ratio
(m (in)) Backcalculation (MPa (ksi))
AC
0.35
0.152 (6.0)
2020 (293)
Base
0.30
0.508 (20.0)
95.8 (13.9)
Subgrade 0.45
Infinity
178.6 (25.9)

E(t) Backcalculation
(MPa (ksi))
Fig. 101
101.4 (14.7)
142.0 (20.6)

Table 44 shows the elastic modulus of base and subgrade layer from VE backcalculation are
very similar as the MODCOMP5 result. It can be noted that MODCOMP5 results are slightly
more consistent since the pavement structure is the same as that from the previous examples in
terms of base and subgrade materials.
163

deflection (mls)

Deflection history simulation (E(t))
40
35
30
25
20
15
10
5
0
-5

sensor 1
sensor 4
Sensor 7

0

5

10

15

20
25
time (ms)

sensor 2
sensor 5
Sensor 8

30

35

sensor 3
sensor 6
sensor 9

40

45

Figure 100: Input and simulated time-histories by backcalculation of E(t) in L11S4 site
(Solid lines represent the input deflection history, symbols represent the VE simulation)

Modulus (Pa)

Relaxation Modulus Comparison: Backcalculated
Vs. Input
1.E+10
Backcalculated
Actual

1.E+09

1.E+08

1.E+07
1.E-03

1.E-02

1.E-01

1.E+00 1.E+01
Time (s)

1.E+02

Figure 101: Input and backcalculated E(t) in L11S4 site

164

1.E+03

Fig. 100 shows that the predicted deflection time histories match very well with the actual
values for all the sensors. Moreover, the backcalculated and input E(t) matches reasonably well
(Fig. 101), especially during the concerned time range 0.0001 – 0.1 s.

6.2.3.4 Discussion
Field FWD tests results indicate the backcalculation algorithm is appropriate in some cases,
especially for the concerned reduced time range 0.0001 – 0.1 s. In other cases, the VE
backcalculation algorithm does not work well, because of the data measurement issues as well as
the insufficient information on the temperature profile of the AC layer. It is noted that the
temperature of the AC layer can have a significant effect on the response and thus the
backcalculated values. If better data acquisition and better interpretation of the FWD test data are
available, the backcalculation results should be better (Kutay et al, 2011).

6.3 Sensitivity Analysis for the AC Layer
The base and subgrade moduli are very sensitive in the backcalculation program, and the
time efficiency is low if the base and subgrade moduli are calculated in the VE backbacalculation
algorithm. Thus, given that the focus of this dissertation is the backcalculation of E(t) and |E*|,
the base and subgrade moduli are assumed to be known, and only AC unknown function is
considered in the sensitivity analysis. The base modulus is typically assumed higher than the
subgrade modulus, to simulate realistic pavement structures.
Four hypothetical pavement structures are given in Table 45. The effect of errors from
deflection measurement and sampling time, and layer thickness are considered for the sensitivity
analysis in the VE backcalculation program. The sampling time is investigated in terms of (1)
165

time shift of the entire deflection impulse, and (2) expansion/contraction of the deflection pulse
in the post-peak zone. Illustrations of both errors in sampling time are shown below in Fig. 102
and Fig. 103. Both dynamic modulus mastercurve |E*| and relaxation modulus E(t) for the AC
layer are independently backcalculated, while assuming the elastic moduli of the base and
subgrade layers are known.

Table 45: Layer information of each pavement structure for sensitivity analysis

Physical
Layer
AC
Base Layer
Subgrade Layer

Elastic Modulus
(MPa (ksi))
Experimental Data
206.9 (30.0)
34.5 (5.0)

Thickness
(m (in))
0.152 (6.0)
0.610 (24.0)
Infinity

B

AC
Base Layer
Subbase Layer
Subgrade Layer

Experimental Data
344.8 (50.0)
137.9 (20.0)
68.9 (10.0)

0.178 (7.0)
0.152 (6.0)
1.143 (45.0)
Infinity

0.35
0.40
0.40
0.45

C

AC
Base Layer
Subgrade Layer
Bedrock

Experimental Data 0.203m (8.0)
344.8 (50.0)
0.381 (15.0)
68.9 (10.0)
1.397 (55.0)
6895.0 (1000.0)
Infinity

0.35
0.40
0.45
0.25

D

AC
CTB Base Layer
Subgrade Layer

Experimental Data
4998.9 (725.0)
68.9 (10.0)

0.35
0.20
0.45

Case
A

166

0.102 (4.0)
0.152 (6.0)
Infinity

Poisson’s
Ratio
0.35
0.30
0.45

Time Offset (1.4ms)

Deflection (inch)

0.03
0.02

actual
Time (+1.4ms)
Time (-1.4ms)

0.02
0.01
0.01
0.00
0.00

0.01

0.02

0.03
0.04
Time (s)

0.05

0.06

0.07

Figure 102: Error caused by time shift of entire pulse for one sensor

Postpeak Time Offset (1.4ms)
Deflection (inch)

0.03
0.02

actual
postpeak (+1.4ms)
postpeak (-1.4ms)

0.02
0.01
0.01
0.00
0.00

0.01

0.02

0.03
0.04
Time (s)

0.05

0.06

0.07

Figure 103: Error caused by expansion/contraction of deflection pulse in the post-peak
zone for one sensor

167

6.3.1 Sensitivity Analysis for Deflection Measurement of the First Sensor
Since the first sensor which is under the FWD test load affects the result of the AC modulus
more, only +2% or -2% error of the deflection history for sensor 1 is considered. The
backcalulation result with inaccurate input is checked with the backcalculation result with
accurate input, as shown in Fig. 104 - 111.

Dynamic Modulus (MPa)

1.E+05

|E*| Mastercurve input vs. backcalculation &
combination

1.E+04
1.E+03

Actual
Back-calculated
+2% measurement
-2% measurement

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Frequency (Hz)
Figure 104: Sensitivity analysis of |E*| for deflection measurement error (Case A)

168

Modulus (Pa)

Relaxation modulus E(t) input vs. backcalculated
1.E+11
Back-Calculated
Actual
1.E+10
+2% measurement
-2% measurement
1.E+09
1.E+08
1.E+07
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 105: Sensitivity analysis of E(t) for deflection measurement error (Case A)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05
1.E+04
1.E+03

Actual
Back-calculated
+2% measurement
-2% measurement

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 106: Sensitivity analysis of |E*| for deflection measurement error (Case B)

169

Modulus (Pa)

Relaxation modulus E(t) input vs. backcalculated
1.E+11
Back-Calculated
Actual
1.E+10
+2% measurement
-2% measurement
1.E+09
1.E+08
1.E+07
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 107: Sensitivity analysis of E(t) for deflection measurement error (Case B)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05
1.E+04
1.E+03

Actual
Back-calculated
+2% measurement
-2% measurement

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 108: Sensitivity analysis of |E*| for deflection measurement error (Case C)

170

Modulus (Pa)

Relaxation modulus E(t) input vs. backcalculated
1.E+11
Back-Calculated
Actual
1.E+10
+2% measurement
-2% measurement
1.E+09
1.E+08
1.E+07
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 109: Sensitivity analysis of E(t) for deflection measurement error (Case C)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05
1.E+04
1.E+03
Actual
Back-calculated
+2% measurement
-2% measurement

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 110: Sensitivity analysis of |E*| for deflection measurement error (Case D)

171

Modulus (Pa)

Relaxation modulus E(t) input vs. backcalculated
1.E+11
Back-Calculated
Actual
1.E+10
+2% measurement
-2% measurement
1.E+09
1.E+08
1.E+07
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 111: Sensitivity analysis of E(t) for deflection measurement error (Case D)

The above cases show the VE backcalculation program is not sensitive to the error of
deflection measurement for sensor 1, especially for the range of interest, i.e., reduced frequency
10 – 10000 Hz in |E*| or reduced time 0.0001 – 0.1 second in E(t). Thus, the measurement of the
deflection in sensor 1 is not significant for the backcalculation result.

6.3.2 Sensitivity Analysis for Time Shift Error of the First Sensor
Still sensor 1 is considered for the sensitivity analysis, for the same reason stated above, with
+1.4 ms or -1.4 ms (depending on sampling time) time shift of the deflection history. The
backcalulation results with inaccurate input are checked with the backcalculation results with
accurate input, as shown in Fig. 112 - 119.
172

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05
1.E+04
1.E+03
Actual
Back-calculated
+1.4ms measurement
-1.4ms measurement

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 112: Sensitivity analysis of |E*| for time shift error of sensor 1 (Case A)

Modulus (Pa)

Relaxation modulus E(t) input vs. backcalculated
1.E+11
1.E+10
1.E+09
1.E+08

Back-Calculated
Actual
+1.4ms measurement
-1.4ms measurement

1.E+07
1.E-05

1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 113: Sensitivity analysis of E(t) for time shift error of sensor 1 (Case A)

173

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05
1.E+04
1.E+03
1.E+02

Actual
Back-calculated
+1.4ms measurement
-1.4ms measurement

1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 114: Sensitivity analysis of |E*| for time shift error of sensor 1 (Case B)

Modulus (Pa)

Relaxation modulus E(t) input vs. backcalculated
1.E+11
1.E+10
1.E+09
1.E+08
1.E+07
1.E-05

Back-Calculated
Actual
+1.4ms measurement
-1.4ms measurement
1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 115: Sensitivity analysis of E(t) for time shift error of sensor 1 (Case B)

174

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05
1.E+04
1.E+03
1.E+02

Actual
Back-calculated
+1.4ms measurement
-1.4ms measurement

1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 116: Sensitivity analysis of |E*| for time shift error of sensor 1 (Case C)

Modulus (Pa)

Relaxation modulus E(t) input vs. backcalculated
1.E+11
1.E+10
1.E+09
1.E+08
1.E+07
1.E-05

Back-Calculated
Actual
+1.4ms measurement
-1.4ms measurement
1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 117: Sensitivity analysis of E(t) for time shift error of sensor 1 (Case C)

175

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05
1.E+04
1.E+03
1.E+02

Actual
Back-calculated
+1.4ms measurement
-1.4ms measurement

1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 118: Sensitivity analysis of |E*| for time shift error of sensor 1 (Case D)

Relaxation modulus E(t) input vs. backcalculated
1.E+12

Modulus (Pa)

1.E+11
1.E+10
1.E+09

Back-Calculated
Actual

1.E+08
1.E+07
1.E-05

+1.4ms measurement
-1.4ms measurement
1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 119: Sensitivity analysis of E(t) for time shift error of sensor 1 (Case D)

176

The above cases show the backcalculation program is significantly sensitive to the time shift
error of sensor 1. Thus, there should not be any error in time recording during the field FWD test.

6.3.3 Sensitivity Analysis for Post-Peak Time Offset of the First Sensor
Still sensor 1 is considered for the sensitivity analysis, for the same reason stated above, with
+1.4 ms or -1.4 ms (depending on sampling time) time shift after the peak deflection occurred.
The backcalulation results with inaccurate input are checked with the backcalculation results
with accurate input, as shown in Fig. 120 - 127.

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05
1.E+04
1.E+03
1.E+02

Actual
Back-calculated
+1.4ms postpeak
-1.4ms postpeak

1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 120: Sensitivity analysis of |E*| for post-peak time offset of sensor 1 (Case A)

177

Modulus (Pa)

Relaxation modulus E(t) input vs. backcalculated
1.E+11
1.E+10
1.E+09
1.E+08
1.E+07
1.E-05

Back-Calculated
Actual
+1.4ms postpeak
-1.4ms postpeak
1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 121: Sensitivity analysis of E(t) for post-peak time offset of sensor 1 (Case A)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05
1.E+04
1.E+03
1.E+02
Actual
Back-calculated
+1.4ms postpeak
-1.4ms postpeak

1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 122: Sensitivity analysis of |E*| for post-peak time offset of sensor 1 (Case B)

178

Modulus (Pa)

1.E+11

Relaxation modulus E(t) input vs. backcalculated

1.E+10
1.E+09
1.E+08
1.E+07
1.E-05

Back-Calculated
Actual
+1.4ms postpeak
-1.4ms postpeak
1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 123: Sensitivity analysis of E(t) for post-peak time offset of sensor 1 (Case B)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
1.E+05
combination
1.E+04
1.E+03
1.E+02
Actual
Back-calculated
+1.4ms postpeak
-1.4ms postpeak

1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 124: Sensitivity analysis of |E*| for post-peak time offset of sensor 1 (Case C)

179

Modulus (Pa)

Relaxation modulus E(t) input vs. backcalculated
1.E+11
1.E+10
1.E+09
1.E+08
1.E+07
1.E-05

Back-Calculated
Actual
+1.4ms postpeak
-1.4ms postpeak
1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 125: Sensitivity analysis of (t) for post-peak time offset of sensor 1 (Case C)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
1.E+05
combination
1.E+04
1.E+03
Actual
Back-calculated
+1.4ms postpeak
-1.4ms postpeak

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 126: Sensitivity analysis of |E*| for post-peak time offset of sensor 1 (Case D)

180

Relaxation modulus E(t) input vs. backcalculated
1.E+11

Modulus (Pa)

1.E+10
1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

Back-Calculated
Actual
+1.4ms postpeak
-1.4ms postpeak
1.E-03

1.E-01
Time (s)

1.E+01

1.E+03

Figure 127: Sensitivity analysis of E(t) for post-peak time offset of sensor 1 (Case D)

Above cases show the backcalculation program is significantly sensitive to the post-peak
time offset of sensor 1, although it works well within the interested range in some cases. Thus,
there cannot be any offset in time recording after the maximum deflection occurs during the field
FWD test.

6.3.4 Sensitivity Analysis for Layer Thickness
Layer thickness is known to be very sensitive in MODCOMP5, and therefore it is expected
here it is as sensitive for this VE backcalculation program. In this analysis, the thickness of AC
layer and base layer were varied by +0.5 in or -0.5 in individually, as well as combined. Each
pavement structure may therefore have 8 cases of layer thickness for a total of 32
backcalculation runs. The backcalulation results with inaccurate input are checked with the
backcalculation results with accurate input, as shown in Fig. 128 - 135.
181

1.E+04
1.E+03
Actual
+0.5 AC
+0.5 BS
+0.5 AC +0.5 BS
-0.5 AC +0.5 BS

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

Back-calculated
-0.5 AC
-0.5 BS
+0.5 AC -0.5 BS
-0.5AC -0.5 BS

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 128: Sensitivity analysis of |E*| for layer thickness (Case A)

1.E+11

Relaxation modulus E(t) input vs.
backcalculated

1.E+10
Modulus (Pa)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05

1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

Back-Calculated
AC+0.5
BS +0.5
AC +0.5, BS +0.5
AC -0.5, BS +0.5
1.E-03

1.E-01
Time (s)

Actual
AC -0.5
BS -0.5
AC +0.5, BS -0.5
AC -0.5, BS -0.5
1.E+01

1.E+03

Figure 129: Sensitivity analysis of E(t) for layer thickness (Case A)

182

1.E+04
1.E+03
Actual
+0.5 AC
+0.5 BS
+0.5 AC +0.5 BS
-0.5 AC +0.5 BS

1.E+02
1.E+01
1.E+00
1.E-06

1.E-04

Back-calculated
-0.5 AC
-0.5 BS
+0.5 AC -0.5 BS
-0.5AC -0.5 BS

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 130: Sensitivity analysis of |E*| for layer thickness (Case B)

1.E+11

Relaxation modulus E(t) input vs.
backcalculated

1.E+10
Modulus (Pa)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
1.E+05
combination

1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

Back-Calculated
AC+0.5
BS +0.5
AC +0.5, BS +0.5
AC -0.5, BS +0.5
1.E-03

1.E-01
Time (s)

Actual
AC -0.5
BS -0.5
AC +0.5, BS -0.5
AC -0.5, BS -0.5
1.E+01

1.E+03

Figure 131: Sensitivity analysis of E(t) for layer thickness (Case B)

183

1.E+04
1.E+03
1.E+02
1.E+01
1.E+00
1.E-06

Actual
+0.5 AC
+0.5 BS
+0.5 AC +0.5 BS
-0.5 AC +0.5 BS
1.E-04

Back-calculated
-0.5 AC
-0.5 BS
+0.5 AC -0.5 BS
-0.5AC -0.5 BS

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 132: Sensitivity analysis of |E*| for layer thickness (Case C)

1.E+11

Relaxation modulus E(t) input vs.
backcalculated

1.E+10
Modulus (Pa)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05

1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

Back-Calculated
AC+0.5
BS +0.5
AC +0.5, BS +0.5
AC -0.5, BS +0.5
1.E-03

Actual
AC -0.5
BS -0.5
AC +0.5, BS -0.5
AC -0.5, BS -0.5

1.E-01
Time (s)

1.E+01

1.E+03

Figure 133: Sensitivity analysis of E(t) for layer thickness (Case C)

184

1.E+04
1.E+03
1.E+02
1.E+01
1.E+00
1.E-06

Actual
+0.5 AC
+0.5 BS
+0.5 AC +0.5 BS
-0.5 AC +0.5 BS
1.E-04

Back-calculated
-0.5 AC
-0.5 BS
+0.5 AC -0.5 BS
-0.5AC -0.5 BS

1.E-02
1.E+00
Frequency (Hz)

1.E+02

1.E+04

Figure 134: Sensitivity analysis of |E*| for layer thickness (Case D)

1.E+11

Relaxation modulus E(t) input vs.
backcalculated

1.E+10
Modulus (Pa)

Dynamic Modulus (MPa)

|E*| Mastercurve input vs. backcalculation &
combination
1.E+05

1.E+09
1.E+08
1.E+07
1.E+06
1.E-05

Back-Calculated
AC+0.5
BS +0.5
AC +0.5, BS +0.5
AC -0.5, BS +0.5
1.E-03

Actual
AC -0.5
BS -0.5
AC +0.5, BS -0.5
AC -0.5, BS -0.5

1.E-01
Time (s)

1.E+01

1.E+03

Figure 135: Sensitivity analysis of E(t) for layer thickness (Case D)

185

Above examples show the backcalculation program can catch the tendency very well, but
the error of the dynamic modulus may vary as high as 100%. Therefore, the thickness of each
physical layer should be as accurate as possible in the field FWD test.

6.3.5 Discussion of Sensitivity Analyses
Above sensitivity analysis show that the time information is extremely important in the VE
backcalculation; great effort is required to make sure the sensors data are sampled at the right
time with the right interval. The accuracy of deflection measurement is important, but not as
significant as time information. The layer thickness can significantly affect the VE
backcalculation results, and it must be very accurate in the field FWD test. The best way of
estimating the layer thickness is field coring or using ground penetration radar (GPR).

6.4 Summary
Four numerical examples from VE forward solution, and four numerical examples from
dynamic solution SAPSI, show the VE backcalculation algorithm works much better
than the static backcalculation result from MODCOMP5. The algorithm is inappropriate
if there is shallow bedrock underneath.
The modulus of base or subgrade layer may vary from the actual result, with error
around 10%. If realistic pavement structure is constructed with the base layer stiffer than
the subgrade layer, like the numerical examples from VE forward solution, the error can
be reduced significantly.
Three field FWD test data were input into the backcalculation program. The program

186

convergence was very good; however, the results may be questionable. Because the data
collection in the field test is very critical, the accuracy of the program needs further
verification.
A sensitivity analysis was done for the main input variables. The analysis showed that
errors in the deflection time-histories are the most significant factor affecting the
backcalculated AC relaxation and complex modulus curves. Also, the layer thicknesses
should be as accurate as possible.

187

Chapter 7. Conclusions and Recommendations
7.1 Summary
In this research, a new method for backcalculating viscoelastic (VE) parameters of a
layered flexible pavement system, based on a layered viscoelastic solution and using FWD
deflection time histories, has been developed. The method allows for theoretically
backcalculating the asphalt concrete complex and relaxation modulus curves, and base and
subgrade layer moduli for a three layer system. A four-parameter sigmoid function is prescribed
to describe the variation of |E*| in frequency domain and of E(t) in time domain, respectively.
The total number of unknowns in the inverse problem is six. In order to effectively use the
layered viscoelastic solution, the time delays for the deflection sensors are eliminated since they
are due to wave propagation. On the other hand, the time delays of the farther sensors are used to
estimate a seed value for the subgrade modulus based on wave propagation theory.
The VE backcalculation algorithm is an extension of the solution used in the MICHBACK
program (Harichandran et al., 1994). It uses the modified Newton’s method to obtain a least
squares solution of an over determined set of equations. In the MICHBACK solution, these sets
are real-valued and correspond to the peak deflection values. In this research, the author uses
deflection time histories, or many deflection basins, since the proposed backcalculation scheme
uses a layered VE solution to predict the time-dependent deflection basins. The singular value
decomposition (SVD) with scaling is used to find the inverse of the gradient matrix. To avoid
local minima, a random function inside the backcalculation algorithm is used to generate a
188

different set of values once a solution is found by the Newton’s method. This is done multiple
times so that a number of different potential solutions are obtained. The best candidate from
these solutions is the set with the minimum error for matching the defection time histories.
Finally, the best candidate set of values is input as seed values in MATLAB’s internal function
“fminsearch” which is then used to select the final solution.
The new program was theoretically verified using synthetic data, and was evaluated using
field FWD data. For the theoretical verifications, time histories of FWD surface deflections
generated from the layered viscoelastic solution and from SAPSI were used to verify the
capabilities of the newly developed VE backcalculation program. Various pavement profiles of
different combinations of layer thicknesses and moduli were analyzed. Some profiles included
cases where there was a shallow bedrock or ground water table. In addition a sensitivity analysis
was conducted to investigate the effects of measurement errors in the deflection time histories,
and layer thicknesses.
To evaluate the applicability of the new VE backcalculation program to interpret field tests,
measured deflection time history data from several FWD tests conducted at the FHWA
Turner-Fairbanks Accelerated Loading Facility were used to backcalculate relaxation modulus
curves for the AC layer. These were then compared with the laboratory derived relaxation
modulus curves obtained from laboratory compacted specimens of the same AC mixtures.

7.2 Conclusions
Based on the theoretical verification analysis, the following conclusions are drawn:
• The backcalculation results show very good agreement with the actual modulus values.

The error in the base and subgrade moduli is generally less than 3%. The AC layer

189

relaxation and complex modulus curves match the actual functions within the range of
0.0001 to 0.1 sec, and 10 to 10,000 Hz, respectively.
• When shallow bedrock exists, the backcalculated and actual deflection time histories and

relaxation and complex modulus curves do not match well. The reason for this discrepancy
is that the stiff layer traps the wave energy within the pavement structure, and thus it
cannot be simulated by the layered VE solution.
The verification of the new VE backcalculation solution using field measured deflection time
histories under FWD testing showed very good agreement between backcalculated and
laboratory derived relaxation modulus curves, especially considering that laboratory and field
conditions are known to be different.
Finally, a sensitivity analysis was done for the main input variables. The analysis showed that
errors in the deflection time-histories are the most significant factor affecting the backcalculated
AC relaxation and complex modulus curves. Also, the layer thicknesses should be as accurate as
possible.

7.3 Recommendations
The proposed VE forward solution is inappropriate to simulate the pavement response under
FWD test if there is shallow bedrock underneath. The reason of this discrepancy is that the stiff
layer traps the wave energy within the pavement structure, and thus it cannot be simulated by the
layered VE solution. Therefore, for these cases, a full dynamic solution which takes into account
inertial and wave propagation effects is required. Further research is necessary to improve the
computational efficiency of the dynamic solution.
The proposed method (k-value) to estimate the elastic modulus of the subgrade layer does
190

not work well if there is shallow bedrock underneath. This is because it is based on the Raleigh
wave solution which is applicable to a deep subgrade (half-space). The Love wave solution is
suitable for a stiff layer condition; however, the solution does not lead to a simple relationship
between wave velocities and elastic modulus of the subgrade.
In order to improve the interpretation of field measured FWD test records, it is
recommended that enhancements in data collection accuracy and FWD technology are pursued.
Examples of such enhancements include longer duration loading, better measurement accuracy
in the unloaded time range, and more accurate temperature measurements in the AC layer. Such
enhancements will improve the accuracy of the AC modulus curves in a wider range of time
and/or frequency.
If the time-temperature shift function is known, FWD test can be done at different
temperatures, and backcalculation is done for each temperature. By combining the results at the
different temperatures, the range of |E*| and E(t) can be extended. However, very low
frequencies and/or very long times, Witczak’s or similar regression equations can help in
estimating the values in these ranges.

191

BIBLIOGRAPHY

192

BIBLIOGRAPHY
W. E. A. Acum, L. Fox. ‘Computation of Load Stresses in a Three-Layer Elastic System’,
Geotechnique, Vol. 2, No. 4, pp 293 - 300, 1951.
R. Al-Khoury, A. Scarps, C. Kasbergen, and J. Blaauwendraad. ‘Spectral Element Technique for
Efficient Parameter Identification of Layered Media. Part I: Forward Calculation’, International
Journal of Solids and Structures, Vol. 38, pp 1605 - 1623, 2001.
W. S. Mogawer, A. J. Austerman, M. E. Kutay, and F. Zhou. ‘Utilizing Linear Viscoelastic
Continuum Damage and Fracture Mechanics Approaches to Evaluate Binder Elasticity as a
Simple Performance Test for Fatigue Cracking of Polymer Modified HMA Mixtures’. Road
Materials and Pavement Design, in press, 2010.
G. Y. Baladi. ‘Integrated Material and Structural Design Method for Flexible Pavements’, Vol. 1:
Technical Report, FHEA, 1988.
Gilbert Y. Baladi, Personal Communications. January, 2010.
D. M. Burmister. ‘The Theory of Stresses and Displacement in Layered Systems and Application
to the Design of Airport Runways’, Proceeding Highway Research Board, Vol. 23, 1943.
A. J. Bush. ‘Computer Program BISDEF’, US Army Corps of Engineer Waterways Experiment
Station, Vicksburg. Miss., 1985.
A. J. Bush, D. R. Alexander. ‘Pavement Evaluation Using Deflection Basin Measurements and
Layered theory’, Transportation Research Record, No. 1022, Transportation Research Board,
National Research Council, Washington D. C., pp 16 - 29, 1985.
K. Chatti, A. Iftikhar, and H. B. Kim. ‘Comparison of Energy Based Fatigue Curves for Asphalt
Mixtures Using Cyclic Indirect Tensile and Flexural Tests’. 2nd International Conference on
Engineering Materials, pp. 271 - 282, 2001.
E. Y. G. Chen, and E. Pan. ‘Responses of a transversely isotropic layered half-space to multiple
horizontal loads.’ ASCE Geotechnical Special Publication: Innovations in the Characterization,
Modeling and Simulation of Pavements and Materials, 2007.
Ewan Yuanguo Chen. ‘Viscoelastic Modeling of Flexible Pavement’, Ph.D. dissertation,
University of Akron, Akron, 2009.
E. Y. G. Chen, E. Pan, and R. Green. ‘Surface loading of a multilayered viscoelastic pavement: A
semi-analytical solution.’ Journal of Engineering Mechanics, ASCE, 2009.
S. S. Chen. ‘The Response of Multilayered Systems to Dynamic Surface Loads’, Ph.D.
dissertation, University of California, Berkeley, 1987.

193

Shengkai Chen, Jianming Ling, Zhigang Luo. ‘Stress-dependent Characteristics and Prediction
Model of the Resilient Modulus of Subgrade Soils’, China Civil Engineering Journal, vol. 40, No.
6, pp. 95 – 99 + 104, June 2007.
K. M. Chua, R. L. Lytton. ‘Load Rating of Light Pavement Structures’, Transportation Research
Record, No. 1043, Transportation Research Board, National Research Council, Washington D. C.,
pp 89 - 102, 1985.
Jo Sias Daniel, Y. Richard Kim, Hyun-Jong Lee. ‘Effects of Aging on Viscoelastic Properties of
Asphalt-Aggregate Mixtures’, Transportation Research Record, No. 1630, Transportation
Research Board, National Research Council, Washington D. C., pp 21 - 27, 1998.
J. F. Doyle. ‘Wave Propagation in Structures: Spectral Analysis Using Fast Discrete Fourier
Transforms’. Springer, New York, 1997.
V. S. Endiran. ‘Nonlinear dynamic computer modeling of flexible pavements’, Ph.D. dissertation,
University of Hong Kong, Hong Kong, 1999.
W. Maurice Ewing, Wenceslas S. Jardetzky, Frank Press, ‘Elastic Waves in Layered Media’,
McGraw-Hill Book Company, Inc. New York, Toronto, London, 1957.
Federal Highway Administration, U.S. Department of Transportation. ‘Resilient Modulus of
Unbound Granular Base/Subbase Materials and Subgrade Soils’, Long-Term Pavement
Performance, Protocol 46, http://www.tfhrc.gov/pavement/ ltpp/pdf/p46.pdf, 1996.
V. Ganji, N. Gucunski, and S. Nazarian. ‘Automated inversion procedure for spectral analysis of
surface waves’, J. Geotech. Engrg., ASCE, v. 124, n. 8, p. 757 – 770, 1998.
R. S. Harichandran, T. Mahmood, A. Raab, and G. Y. Baladi. ‘Backcalculation of Pavement
Layer Moduli, Thicknesses and Bedrock Depth Using a Modified Newton Method’, In
Nondestructive Testing of Pavements and Backcalculation of Moduli (Second Volume), ASTM
STP 1198, H. L. Von Quintus, A. J. Bush and G. Y. Baladi (eds.) American Society for Testing
and Materials, Philadelphia, PA, pp. 68 - 82, 1994.
Ronald S. Harichandran, Personal Communications. May, 2008.
H. Y. Huang. ‘Pavement Analysis and Design’, Prentice Hall, Pearson Educaion, Inc., Upper
Saddle River, New Jersey, 2004.
L. H. Irwin. ‘Instructional Guide for Backcalculation and the Use of MODCOMP’, CLRP
Publication No. 94-10, Cornell University, Local Roads Program, Ithaca, NY, March, 1994.
Yigong Ji. ‘Frequency and Time Domain Backcalculation of Flexible Pavement Layer
Parameters’. Ph.D. dissertation, Michigan State University, East Lansing, 2005.
E. Kausel and J. M. Roesset. ‘Stiffness Matrices for Layered Soils’. Bulletin of the Seismological
Society of America, Vol. 71 No. 6, pp.1743 ~ 1761, December, 1981.

194

E. Kausel, R. Peek. ‘Dynamic Loads in the Interior of a Layered Stratum: An Explicit Solution’,
Bulletin of the Seismological Society of America, Vol. 72, No. 5, pp 1459 - 1482, 1982.
Y. Richard Kim, Youngguk Seo, Mark King, and Mostafa Momen. ‘Dynamic Modulus Testing of
Asphalt Concrete in Indirect Tensile Mode’. Transportation Research Record, n 1891, pp. 163 173, 2004.
Y. Richard Kim. ‘Modeling of Asphalt Concrete’, McGraw-Hill, August, 2008.
Y. Richard Kim, N. Guddati Murthy, B. Shane Underwood, et al. ‘Development of a Multiaxial
VEPCD-FEP++’. Final report, May, 2008.
M. E. Kutay, N. Gibson and J. Youtchef. On-going unpublished research involving analysis of
FHWA ALF data, 2008 - 2009.
M. E. Kutay, K. Chatti, L. Lei. ‘Backcalculation of Dynamic Modulus (|E*|) Mastercurve from
FWD Surface Deflections’. Transportation Research Board, Proceeding, Washington D.C., 2011.
E. Levenberg, ‘Validation of NCAT Structural Test Track Experiment using INDOT APT Facility’,
Final Report. North Central Superpave Center, Joint Transportation Research Program Project
No. C-36-31R SPR-2813, Purdue University, 2008.
T. Mahmood. ‘Backcalculation of Pavement Layer Properties from Deflection Data’, Ph.D.
dissertation, Michigan State University, East Lansing, 1993.
NCHRP. ‘2002 Design Guide: Design of New and Rehabilitated Pavement Structures’, Draft
Final Report, NCHRP Study 1-37A, National Cooperative Highway Research Program,
Washington D.C., July, 2004.
A. T. Papagiannakis and E.A. Masad. ‘Pavement Design and Materials’, John Wiley & Sons,
Inc., Hoboken, New Jersey, 2008.
S.W. Park and R.A. Schapery, ‘Methods of Interconversion between Linear Viscoelastic Material
Functions. Part I – a Numerical Method Based on Prony Series.’, International Journal of Solids
and Structures, Vol. 36, pp. 1653 - 1675, 1999.
W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling. ‘Numerical Recipes’, Cambridge
University Press, Cambridge, 1989.
Lufti Raad, Jose L. Figueroa. ‘Load Response of Transportation Support System’, Journal of
Transportation Engineering, American Society of Civil Engineers, Vol. 106, No. 1, pp 111 - 128,
1980.
F. E. Richart, Jr., R. D. Woods, J. R. Hall, Jr. ‘Vibrations of Soils and Foundations’, Prentice-Hall,
Inc., Englewood Cliffs, New Jersey. 1970.

195

S. A. Rizzi, J. F. Doyle. ‘A spectral element approach to wave motion in layered solids’. Journal
of Vibration, Acoustics, Stress, and Reliability in Design. Vol. 114, No. 2, pp 133 - 140, April,
1992.
J. M. Roesset, Kenneth H. Stokoe II, Chia-Ray Seng. ‘Determination of Depth to Bedrock from
Falling Weight Deflectometer Test Data’, Transportation Research Record, No. 1504,
Transportation Research Board, National Research Council, Washington D. C., pp 68 - 78, 1995.
R.A. Schapery. ‘Method of Viscoelastic Stress Analysis using Elastic Solutions’. Journal of the
Franklin Institute, Vol. 279, No. 4, pp 268 ~ 289, 1965.
R.A. Schapery. ‘Viscoelastic Behavior and Analysis of Composite Materials’. Mechanics of
Composite Materials. New York, Academic Press, Inc., pp 85 - 168, 1974.
R. A. Schapery. ‘Correspondence Principles and a Generalized J Integral for Large Deformation
and Fracture Analysis of Viscoelastic Media’. International Journal of Fracture, No. 25, pp 195 223, 1984.
Richard N. Stubstad, Shiraz D. Tayabji, Erland O. Lukanen. ‘LTPP Data Analysis: Variations in
Pavement Design Inputs’, NCHRP Web Document 48 (Project 20-50[5]),
http://onlinepubs.trb.org/onlinepubs/nchrp/nchrp_w48.pdf, 2002.
Richard Stubstad, Jiang Y. Jane, Erland Lukanen. ‘Forwardcalculation of Pavement Moduli with
Load-Deflection Data’, Transportation Research Record, No. 2005, pp 104 - 111, 2007.
J. Uzan. ‘Dynamic Linear Backcalculation of Pavement Material Parameters’, Journal of
Transportation Engineering, American Society of Civil Engineers, Vol. 120, No. 1, pp 109 - 126,
1994.
H. Warren, W. L. Dieckmann. ‘Numerical Computation of Stresses and Strains in a
Multiple-layer Asphalt Pavement System’, Internal report, Chevron Research Corporation,
Richmond, AC, 1963.
Ning Xiang, Personal Communications. November, 2009.
M. S. Yeh. ‘Nonlinear Finite Element Analysis and Design of Flexible Pavements’. Ph.D.
Dissertation. Michigan State University, East Lansing, 1989.
E. J. Yoder, M. W. Witczak, ‘Principles of Pavement Design’, Wiley, New York, 1975.

196