FREEZING AND THAWING OF FROST-SUSCEPTIBLE SOILS (DEVELOPMENT OF A RELIABLE PREDICTIVE MODEL) By Pegah Rajaei A DISSERTATION Submitted to Michigan State University in partial fulfillment of requirements for the degree of Civil Engineering- Doctor of Philosophy 2017 ABSTRACT FREEZING AND THAWING OF FROST-SUSCEPTIBLE SOILS (DEVELOPMENT OF A RELIABLE PREDICTIVE MODEL) By Pegah Rajaei Frost depth is an important factor that affects the design of various transportation infrastructures including pavements, retaining structures, bridge foundations, utility lines, and so forth. Soil freezing can lead to frost heave and heave pressure, which may cause serious stability issues. On the other hand, at the beginning of spring season, the ice starts to thaw from the top down and to a lesser extend from the bottom up. The melted water below the pavement surface is trapped (setting on impermeable frozen materials). It saturates the top part of the upper pavement layer. Consequently, the stiffness of the saturated layer decreases causing substantial decrease in its load bearing capacity and high deformations, which lead to premature and localized failure. To decrease the spring thaw damage, Spring Load Restrictions (SLR) signs are usually placed along the roads. The objectives of this study are to develop accurate and reliable frost and thaw depth and frost heave prediction models, estimate heave pressure and develop a reliable SLR policy. After extensive literature review, various existing frost depth models were identified and tested. These include the finite difference UNSAT-H, the Stefan, the Modified Berggren, and the Chisholm and Phang models. Unfortunately, some of these models require substantial input data that are not available and all models yielded inaccurate results. Therefore, statistical frost depth models were developed using frost depth and air temperature data collected by Michigan Department of Transportation (MDOT); one for clayey soils and one for sandy soil. The two models were then combined using the measured thermal conductivity of clayey and sandy soils. The combined statistical model was then verified using frost depth and air temperature data collected by Minnesota Department of Transportation (MnDOT). Additionally, The Gilpin’s mechanistic-empirical model was employed to predict frost heave. The model produced inaccurate and counterintuitive results in some cases. Therefore, the model was modified and the empirical frost depth model developed in this study was incorporated into the model. The resulting model was then simplified to replace some of the required of input data that are not available. The modified model accuracy was assessed using the frost heave data measured at 5 sites in Oakland County, Michigan. Further, the relationship between frost heave and heave pressures were established for four soil types. Moreover, a new statistical model was developed for calculating the cumulative thaw degree-day (CTDD) using pavement surface temperature and air temperate data collected by MDOT. Then, the thaw depth data measured in the state of Michigan were used to assess Nixon and McRoberts thaw depth predictions model. Since the model did not produce accurate and acceptable results, statistical thaw depth models were developed using the calculated CTDD values and thaw depth data collected by MDOT and MnDOT; one for clayey soils and one for sandy soils. The models were then verified using the calculated CTDD values and thaw depth data collected by MnDOT. Finally, based on the results of thaw depth model a new SLR policy was proposed. In dedication to my beloved parents, Without whom none of my success would be possible. iv ACKNOWLEDGEMENTS I would like to express my special appreciation and thanks to my advisor, Dr. Gilbert Baladi, for providing me the opportunity to complete this work, encouraging my research and for allowing me to grow as a research scientist. His advice on both research as well as on my career have been priceless. I would also like to thank my committee members, Dr. Karim Chatti, Dr. Emin Kutay and Dr. Neil Wright for serving as members of my PhD committee, and for their professional comments and suggestions. I would also like to thank the Michigan Department of Transportation (MDOT) for the financial support and the MDOT staff for their cooperation and for providing soil samples and the soil temperature data. In addition, I gratefully acknowledge the Minnesota Department of Transportation (MNDOT) for providing the soil temperature data and for its assistance. At the end, a special thanks to my family. Words cannot express how grateful I am to my sister, brother, mother, and father for all of the sacrifices that they have made on my behalf. I would also like to thank all of my friends, particularly Nan Hu, Andisheh Ranjbari and Jenifer Banks who supported me in writing, and incented me to strive towards my goal. v TABLE OF CONTENTS LIST OF TABLES ....................................................................................................................... viii LIST OF FIGURES ...................................................................................................................... xii CHAPTER 1 ................................................................................................................................... 1 INTRODUCTION & RESEARCH PLAN .................................................................................... 1 1.1 Problem Statement ........................................................................................................... 1 1.2 Study Objectives .............................................................................................................. 3 1.3 Research Plan ................................................................................................................... 3 1.3.1 Task 1 - Conduct Comprehensive Literature Review ............................................... 3 1.3.2 Task 2 - Development of Heat Transfer Predictive Model ....................................... 4 1.3.3 Task 3 – Development of Coupled Heat and Mass Transfer Models for Prediction of Frost Heave and Frost Pressure ........................................................................................... 4 1.3.4 Task 4 – Development of SLR Policy ...................................................................... 5 1.4 Dissertation Layout .......................................................................................................... 5 CHAPTER 2 ................................................................................................................................... 7 LITERATURE REVIEW ............................................................................................................... 7 2.1 Frost Depth ....................................................................................................................... 7 Numerical Models ..................................................................................................... 7 2.2.2.1 UNSAT-H Modeling ........................................................................................... 14 Mechanistic Empirical Models ............................................................................... 15 Empirical Models .................................................................................................... 21 2.2 Frost Heave .................................................................................................................... 24 Frost Heave Mitigation ........................................................................................... 24 2.2.2.1 Capillary Theory ................................................................................................. 28 2.2.2.1 Frozen Fringe Theory .......................................................................................... 32 Frost Pressure .......................................................................................................... 51 2.3 Thaw Depth and Seasonal Load Restrictions ................................................................. 54 Thaw Depth ............................................................................................................. 54 Seasonal Load Restriction....................................................................................... 56 CHAPTER 3 ................................................................................................................................. 72 DATA MINING ........................................................................................................................... 72 3.1 Data Base........................................................................................................................ 72 3.2 Frost Depth Data ............................................................................................................ 73 3.2.1 The State of Michigan ............................................................................................. 73 3.2.2 The State of Minnesota ........................................................................................... 79 3.3 Soil Properties Data ........................................................................................................ 79 3.4 Frost Heave Data ............................................................................................................ 82 3.5 Pavement Surface Temperature and Thaw Depth Data ................................................. 83 vi CHAPTER 4 ................................................................................................................................. 85 DATA ANALYSIS & DISCUSSION .......................................................................................... 85 Introduction and Research Objectives............................................................................ 85 Hypotheses ..................................................................................................................... 85 Frost Depth ..................................................................................................................... 87 4.3.1 UNSAT-H Model.................................................................................................... 88 4.3.2 Freezing Index and Freezing Degree Day Calculation ........................................... 90 4.3.2.1 Minnesota Cumulative Freezing Degree Day ..................................................... 90 4.3.2.2 Boyd Cumulative Freezing Degree Day ............................................................. 91 4.3.3 Existing Frost Depth Prediction Models ................................................................. 94 4.3.3.1 Stefan’s Equation ................................................................................................ 94 4.3.3.2 Modified Berggren’s Equation ............................................................................ 97 4.3.3.3 Chisholm’ and Phang’s Equation ...................................................................... 100 4.3.4 Frost Depth Empirical Models .............................................................................. 104 Frost Heave .................................................................................................................. 121 4.4.1 Insulation Effect on the Frost Depth ..................................................................... 126 4.4.1.1 Example............................................................................................................. 129 4.4.2 Gilpin’s Frost Heave Model.................................................................................. 130 4.4.2.1 Basic Assumptions ............................................................................................ 130 4.4.2.2 Heat and Mass Balance Equations .................................................................... 132 4.4.2.3 Ice Pressure Distribution in the Frozen Fringe Zone ........................................ 134 4.4.3 Revised Frost Heave Model .................................................................................. 136 4.4.4 Discussion of the Results of the Revised Frost Heave Model .............................. 138 4.4.5 Heave Pressure ...................................................................................................... 146 Thaw Depth .................................................................................................................. 150 4.5.1 Calculation of Cumulative Thawing Degree day (CTDD) ................................... 150 4.5.1.1 Pavement Surface Modeling ............................................................................. 151 4.5.1.2 Model Verification ............................................................................................ 155 4.5.2 Nixon and McRoberts Equation............................................................................ 160 4.5.3 Thaw Depth Empirical Models ............................................................................. 161 Spring Load Restrictions .............................................................................................. 169 CHAPTER 5 ............................................................................................................................... 172 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS .............................................. 172 5.1 Summary ...................................................................................................................... 172 5.1.1 Frost Depth Modeling ........................................................................................... 172 5.2 Frost Heave Model ....................................................................................................... 173 5.2.1 Thaw Depth Model ............................................................................................... 173 5.3 Conclusions .................................................................................................................. 174 5.4 Recommendation .......................................................................................................... 176 APPENDICES………………………………...………………………………………………..177 APPENDIX A: Frost And Thaw Depth Analysis ……………………………………....…178 APPENDIX B: Frost Heave Stations Profile………………………………………………....198 REFERENCES …………………………………………………………………………….…..205 vii LIST OF TABLES Table 2.1 Thermal resistance values (R-values) at different mean temperature........................... 26 Table 2.2 Design values for FPSF insulation materials based on ACSE 32-01 ........................... 27 Table 2.3 Yearly variation in the resilient modulus of roadbed soils and the associated damage (AASHTO, 1993) .......................................................................................................................... 58 Table 2.4 CTDD threshold for posting SLR (Mahoney et.al., 1986) ........................................... 60 Table 2.5 Complementary guild for SLR .................................................................................... 62 Table 2.6 Reference temperature .................................................................................................. 64 Table 3.1 RWIS stations in the State of Michigan........................................................................ 75 Table 3.2 Data availability in each database ................................................................................. 78 Table 3.3 Sensor use guide (Decagon Devices Inc., 2008). ......................................................... 81 Table 3.4 Measured thermal properties of different types of soil using KD2 Pro ........................ 83 Table 3.5 Measured total heave and frost depths in different soil types, I75, Oakland County, Michigan. ...................................................................................................................................... 84 Table 4.1 Hydraulic and thermal properties for different layers in UNSAT-H model ................. 89 Table 4.2 Cumulative freezing degree day calculation, Waters station, Lower Peninsula ........... 91 Table 4.3 Maximum frost depth predicted by Stefan’s equation for RWIS stations .................... 98 Table 4.4 Frost depth calculation using the modified Berggren’s equation, Benzonia, Lower Peninsula, Michigan .................................................................................................................... 101 Table 4.5 Maximum frost depth predicted by Modified Berggren’s equation for RWIS stations ..................................................................................................................................................... 103 Table 4.6 Performance metrics of Equation 4.8 and Equation 4.9 for frost depth estimations in the State of Michigan and Minnesota ......................................................................................... 108 Table 4.7 Performance metrics of Equation 4.8 and Equation 4.9 for frost depth estimations in the State of Michigan and Minnesota ......................................................................................... 115 Table 4.8 Performance metrics of Equation 4.12 for frost depth estimations in the State of Michigan and Minnesota............................................................................................................. 119 viii Table 4.9 Frost susceptibility classification (COE, 1984) .......................................................... 124 Table 4.10 Different input values for each site, I75, Oakland County, Michigan. ..................... 145 Table 4.11 Different input values for each soil type ................................................................... 147 Table 4.12 Statistical coefficients in Equation 4-30 for each soil type. ..................................... 150 Table 4.13 Performance metrics for the model ........................................................................... 153 Table 4.14 Performance metrics of Equation 4.42, MNDOT and WSDOT models for the average daily pavement surface temperature for 12 RWIS stations in the LP and 9 RWIS stations in the UP in 2012, 2014, and 2015........................................................................................................ 159 Table 4.15 Performance metrics of Equation 4.44, MnDOT and WSDOT models for the CTDD calculation for 12 RWIS stations in the LP and 9 RWIS stations in the UP in 2012, 2014, and 2015............................................................................................................................................. 159 Table 4.16 Maximum thaw depth predicted by Nixon and McRoberts equation for RWIS stations ..................................................................................................................................................... 161 Table 4.17 Performance metrics of Equation 4.45 for thaw depth estimations in the State of Michigan and Minnesota............................................................................................................. 165 Table 4.18 Performance metrics of Equation 4.46 for thaw depth estimations in the State of Minnesota.................................................................................................................................... 169 Table 4.19 Time gap between the measured and calculated thaw depth for sandy and clayey soil in the state of Minnesota……………………………………………………………………..…170 Table A.1 Frost depth calculation using Stefan equation, Benzonia, LP ....................................180 Table A.2 Frost depth calculation using Modified Berggren equation, Benzonia, LP ................180 Table A.3 Thaw depth calculation using Nixon and McRoberts equation, Benzonia, LP ..........180 Table A.4 Frost depth calculation using Stefan equation, Cadillac, LP ......................................181 Table A.5 Frost depth calculation using Modified Berggren equation, Cadillac, LP ..................181 Table A.6 Thaw depth calculation using Nixon and McRoberts equation, Cadillac, LP ............181 Table A.7 Frost depth calculation using Stefan equation, Grayling, LP .....................................182 Table A.8 Frost depth calculation using Modified Berggren equation, Grayling, LP .................182 Table A.9 Thaw depth calculation using Nixon and McRoberts equation, Grayling, LP ...........182 ix Table A.10 Frost depth calculation using Stefan equation, Houghton Lake, LP.........................183 Table A.11 Frost depth calculation using Modified Berggren equation, Houghton Lake, LP ....183 Table A.12 Thaw depth calculation using Nixon and McRoberts equation, Houghton Lake, LP183 Table A.13 Frost depth calculation using Stefan equation, Ludington, LP .................................184 Table A.14 Frost depth calculation using Modified Berggren equation, Ludington, LP ............184 Table A.15 Frost depth calculation using Stefan equation, Reed City, LP .................................185 Table A.16 Frost depth calculation using Modified Berggren equation, Reed City, LP .............185 Table A.17 Thaw depth calculation using Nixon and McRoberts equation, Reed City, LP .......185 Table A.18 Frost depth calculation using Stefan equation, Waters, LP ......................................186 Table A.19 Frost depth calculation using Modified Berggren equation, Waters, LP ..................186 Table A.20 Thaw depth calculation using Nixon and McRoberts equation, Waters, LP ............186 Table A.21 Frost depth calculation using Stefan equation, Williamsburg, LP............................187 Table A.22 Frost depth calculation using Modified Berggren equation, Williamsburg, LP .......187 Table A.23 Frost depth calculation using Stefan equation, Au Train, UP...................................188 Table A.24 Frost depth calculation using Modified Berggren equation, Au Train, UP ..............188 Table A.25 Thaw depth calculation using Nixon and McRoberts equation, Au Train, UP ........188 Table A.26 Frost depth calculation using Stefan equation, Brevort, UP .....................................189 Table A.27 Frost depth calculation using Modified Berggren equation, Brevort, UP ................189 Table A.28 Thaw depth calculation using Nixon and McRoberts equation, Brevort, UP ...........189 Table A.29 Frost depth calculation using Stefan equation, Cooks, UP .......................................190 Table A.30 Frost depth calculation using Modified Berggren equation, Cooks, UP ..................190 Table A.31 Frost depth calculation using Stefan equation, Engadine, UP ..................................191 Table A.32 Frost depth calculation using Modified Berggren equation, Engadine, UP .............191 x Table A.33 Frost depth calculation using Stefan equation, Golden Lake, UP ............................192 Table A.34 Frost depth calculation using Modified Berggren equation, Golden Lake, UP ........192 Table A.35 Thaw depth calculation using Nixon and McRoberts equation, Golden Lake, UP ..192 Table A.36 Frost depth calculation using Stefan equation, Harvey, UP .....................................193 Table A.37 Frost depth calculation using Modified Berggren equation, Harvey, UP .................193 Table A.38 Thaw depth calculation using Nixon and McRoberts equation, Harvey, UP ...........193 Table A.39 Frost depth calculation using Stefan equation, Michigamme, UP ............................194 Table A.40 Frost depth calculation using Modified Berggren equation, Michigamme, UP .......194 Table A.41 Frost depth calculation using Stefan equation, Seney, UP ......................................195 Table A.42 Frost depth calculation using Modified Berggren equation, Seney, UP ...................195 Table A.43 Thaw depth calculation using Nixon and McRoberts equation, Seney, UP .............195 Table A.44 Frost depth calculation using Stefan equation, St. Ignace, UP .................................196 Table A.45 Frost depth calculation using Modified Berggren equation, St. Ignace, UP .............196 Table A.46 Frost depth calculation using Stefan equation, Twin Lakes, UP .............................197 Table A.47 Frost depth calculation using Modified Berggren equation, Twin Lakes, UP ..........197 Table A.48 Thaw depth calculation using Nixon and McRoberts equation, Twin Lakes, UP ....197 xi LIST OF FIGURES Figure 2.1 Heat transfer between pavement surface and air temperature (After Dempsey and Pur, 1990) ............................................................................................................................................. 11 Figure 2.2 A schematic representation of two phase heat conduction .......................................... 16 Figure 2.3 Fusion parameter (μ) versus correction factor (λ) ....................................................... 20 Figure 2.4 a. Equilibrium between ice and water, b. simple ice-water model .............................. 30 Figure 2.5 Schematic diagrams for the frost heave process (Peppin and Style, 2012) ................. 32 Figure 2.6 Schematic diagram of a freezing soil with frozen fringe (After Peppin and Style, 2012) ............................................................................................................................................. 33 Figure 2.7 A schematic representations of equilibrium conditions for ice and water near a substrate (Gilpin, 1980) ................................................................................................................ 36 Figure 2.8 the frost heave simulation model (Gilpin, 1980) ......................................................... 38 Figure 2.9 a zone of frozen soil, a freezing fringe, and an underlying zone of unfrozen soil (Nixon ,1991) ................................................................................................................................ 41 Figure 2.10 Numerical verification of linear temperature profile assumption (Nixon 1991) ...... 43 Figure 2.11 Predicted and observed heave for Konrad test No. 4 (Nixon, 1991) ......................... 43 Figure 2.12 Frost action under foundation causing uplift pressure and behind retaining structure causing horizontal pressure ........................................................................................................... 53 Figure 2.13 Frost heave and thaw consolidation process (Dore, 2004) ........................................ 57 Figure 2.14 MnDOT CDTT calculation flowchart. ..................................................................... 63 Figure 2.15 Sinusoidal temperature variation (After Berg et al) .................................................. 66 Figure 3.1 RWIS station locations, Michigan ............................................................................... 74 Figure 3.2 A soil profile at a typical RWIS station showing the depths of the temperature sensors. ....................................................................................................................................................... 76 Figure 3.3 Maximum frost depth contours in a typical year in the State of Michigan. ................ 77 Figure 3.4 MNDOT stations location, Minnesota. ....................................................................... 79 Figure 3.5 Thermal conductivity measurement using KD2 pro. .................................................. 81 xii Figure 3.6 Five locations where the soil thermal properties were measured in each soil sample using KD2 pro. .............................................................................................................................. 82 Figure 3.7 MDOT frost heave station locations, Oakland County, Michigan. ............................. 84 Figure 4.1 Flow chart of the research plan ................................................................................... 86 Figure 4.2 Schematic of the cross section of the two modeled pavement sites ............................ 88 Figure 4.3 Calculated frost depth using UNSAT-H model and measured frost depth versus time in Waters station, Lower Peninsula, Michigan. ............................................................................ 89 Figure 4.4 Calculated frost depth using UNSAT-H model and measured frost depth versus time in Cadillac station, Lower Peninsula, Michigan. .......................................................................... 90 Figure 4.5 Calculation of freezing index using cumulative freezing degree day ......................... 92 Figure 4.6 Comparison of calculated CFDD using Boyd (Boyd 1976) and Minnesota (MnDOT 2009) methods. .............................................................................................................................. 94 Figure 4.7 soil Thermal conductivity of different types of soil based on water content and dry density obtained by US army cold region and engineering laboratory (CRREL) (Edgar, 2014) . 96 Figure 4.8 The maximum frost depths predicted by Stefan equation versus the measured maximum frost depths in Michigan. ............................................................................................. 97 Figure 4.9 Measured maximum frost depths in Michigan versus the maximum calculated ones using the modified Berggren’s equation ..................................................................................... 102 Figure 4.10 Measured maximum frost depths versus calculated ones using Chisholm’ and Phang’s equation. ........................................................................................................................ 104 Figure 4.11 Frost depths versus cumulative freezing degree day for clayey and sandy soils in the State of Michigan. ....................................................................................................................... 106 Figure 4.12 Measured frost depths in State of Michigan versus the calculated ones using Equation 4.8 ................................................................................................................................ 107 Figure 4.13 Measured frost depth in State of Minnesota versus the calculated ones using Equation 4.8 ................................................................................................................................ 107 Figure 4.14 Measured frost depth in State of Michigan versus the calculated ones using Equation 4.9................................................................................................................................................ 109 Figure 4.15 Measured frost depth in State of Minnesota versus the calculated ones using Equation 4.9 ................................................................................................................................ 109 Figure 4.16 Measured frost depths in Michigan versus cumulative freezing degree day for clayey soil showing the best fit and the U.S. Corps of Engineers equations. ........................................ 111 xiii Figure 4.17 Measured versus calculated frost depth data in clayey soils in Michigan (Equation 4.10). ........................................................................................................................................... 111 Figure 4.18 Frost depths versus cumulative freezing degree day for sandy soil showing the best fit and the U.S. Corps of Engineers equations in the State of Michigan. ................................... 112 Figure 4.19 Measured versus calculated frost depths in sandy soils in Michigan (Equation 4.11). ..................................................................................................................................................... 112 Figure 4.20 Measured frost depths in clayey soil in the state of Minnesota versus the frost depth values calculated using Equation 4.10. ....................................................................................... 114 Figure 4.21 Measured frost depths in sandy soil in the state of Minnesota versus the frost depth values calculated using Equation 4.10. ....................................................................................... 114 Figure 4.22 Correlation between the statistical power coefficient (b) of Equation 4.10 and Equation 4.11 and the corresponding average thermal conductivity of the soil ......................... 117 Figure 4.23 Correlation between the statistical coefficient (a) of Equation 4.10 and Equation 4.11 and the corresponding average thermal conductivity of the soil ................................................ 117 Figure 4.24 Calculated frost depths using Equation 4.12 versus the measured frost depth in clayey and sandy soil in the State of Michigan........................................................................... 118 Figure 4.25 Calculated frost depths using Equation 4.12 versus the measured frost depth in clayey and sandy soil in the State of Minnesota ......................................................................... 119 Figure 4.26 Frost depths versus cumulative freezing degree-day for clayey soil showing the best fit statistical model and the proposed model (Equation 4.12) in Minnesota. ............................ 120 Figure 4.27 Frost depths versus cumulative freezing degree-day for sandy soil showing the best fit and proposed model (Equation 4.12) in the State of Minnesota ........................................... 121 Figure 4.28 Effect of capillary and permeability on frost susceptibility (ACPA, 2008). ........... 122 Figure 4.29 Heaving Rate in laboratory test on different disturbed soil types (COE, 1984)...... 123 Figure 4.30 Frost susceptibility criteria, Canadian Department of Transportation (Edgar, 2014). ..................................................................................................................................................... 125 Figure 4. 31 AASHTO four environmental regions (ACPA, 2008) ........................................... 127 Figure 4.32 The schematic of frost heave model (Gilpin, 1980) ................................................ 131 Figure 4.33 Ice pressure along frozen fringe zone (Gilpin, 1980) .............................................. 135 Figure 4.34 Particle separation pressure ..................................................................................... 136 xiv Figure 4.35 Calculated frost heave for three soil types when GWL= 9 m, TTOP= -3 oC for 100 days, POB = 150 kPa. ................................................................................................................... 139 Figure 4.36 Calculated frozen fringe for three soil types when GWL= 9 m, TTOP= -3 oC for 100 days, POB = 150 kPa. ................................................................................................................... 140 Figure 4.37 Calculated total heave versus overburden pressure for clayey silt in different ground water table depths when TTOP= -3 oC in 100 days. ..................................................................... 141 Figure 4.38 Calculated total heave versus overburden pressure for sandy clayey silt in different ground water table depths when TTOP= -3 oC in 100 days. ......................................................... 141 Figure 4.39 Calculated total heave versus overburden pressure for fine sand and silt with pebbles in different ground water table depths when TTOP= -3 oC in 100 days. ...................................... 142 Figure 4.40 Calculated total heave versus overburden pressure for clayey, silty, gravely, sand in different ground water table depths when TTOP= -3 oC in 100 days. .......................................... 142 Figure 4.41 Total heave versus overburden pressure for clayey silt when TTOP is fixed at TTOP= -3 o C and when TTOP is decreasing with a rate of -.057 per day in 100 days. ................................. 143 Figure 4.42 Measured versus calculated frost heave under the shoulder and pavement in 5 sites, Oakland County, Michigan. ........................................................................................................ 145 Figure 4.43 heave pressure versus calculated total heave for clayey silt in different ground water table depths when TTOP= -3 oC in 100 days. ............................................................................... 148 Figure 4.44 Heave pressure versus calculated total heave in four soil types when TTOP= -3 oC over 100 days. ............................................................................................................................. 149 Figure 4.45 Calculated versus measured surface temperature in 2011 using Equation 4.41 and Equation 4.42 for 12 RWIS stations in the first 120 days of the year. The data are shown in 10day interval.................................................................................................................................. 154 Figure 4.46 Calculated versus measured surface temperature in 2012, 2014 and 2015 for 12 RWIS stations in the LP of the State of Michigan. ..................................................................... 156 Figure 4.47 Calculated versus measured surface temperature in 2012, 2014 and 2015 for 9 RWIS stations in the UP of the State of Michigan. .................................................................... 156 Figure 4.48 CTDD of the calculated surface temperature versus CTDD of the measured surface temperature in 2012, 2014, and 2015 for 12 RWIS stations in the LP of the State of Michigan. ..................................................................................................................................................... 158 Figure 4.49 CTDD of the calculated surface temperature versus CTDD of the measured surface temperature in 2012, 2014, and 2015 for 9 RWIS stations in the UP of the State of Michigan. 158 xv Figure 4.50 Maximum thaw depths predicted by Nixon and McRoberts equation versus the measured maximum thaw depths in Michigan. .......................................................................... 160 Figure 4.51 Measured thaw depths versus cumulative thawing degree-day for sandy soils in the State of Michigan ........................................................................................................................ 163 Figure 4.52 Measured thaw depth in sandy soils in the state of Michigan versus the calculated values using Equation 4.45 ......................................................................................................... 163 Figure 4.53 Measured thaw depths in sandy soil in the state of Minnesota versus the thaw depth values calculated using Equation 4.45. ....................................................................................... 164 Figure 4.54 thaw depths versus cumulative thawing degree-day for sandy soil showing the best fit and proposed model (Equation 4.45) in the State of Minnesota. ........................................... 165 Figure 4.55 Thaw depths versus cumulative thawing degree-day for clayey soils in the state of Minnesota.................................................................................................................................... 167 Figure 4.56 Measured thaw depths in clayey soil in the state of Minnesota (Ada, Marshal, Starbucks, and Rochester) versus the thaw depth values calculated using Equation 4.46. ........ 167 Figure 4.57 Measured thaw depths in clayey soil in Gatzke station versus the thaw depth values calculated using Equation 4.46. .................................................................................................. 168 Figure B.1 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 528+88. ..................................................................................................................................................... 200 Figure B.2 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 652+00. ..................................................................................................................................................... 201 Figure B.3 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 719+00. ..................................................................................................................................................... 202 Figure B.4 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 724+00. ..................................................................................................................................................... 203 Figure B.5 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 474+00. ..................................................................................................................................................... 204 xvi CHAPTER 1 INTRODUCTION & RESEARCH PLAN 1.1 Problem Statement In cold regions such as Michigan, where air temperature drops below 0 C for extended periods of time, frost depth is an important factor that affects the design of infrastructures including pavements, building and bridge foundations and/or utility lines. As the winter begins the pavement starts to freeze from the top down. empirical and localized models for predicting frost depth in soils are fairly well established and not complex for a single layered system. However, it is a challenging task for multi-layered systems subjected to various surface boundary conditions where the air temperature fluctuates. Perhaps complex and sophisticated models to predict the propagation of freezing and thawing fronts can be developed. However, such models require a large amount of data that are expensive to obtain and hence, they cannot be easily implemented. Therefore many State Highway Agencies (SHAs) tend to use simplified analytical or empirical models for estimating frost and thaw depth. The problem with such models is that they require calibration from one state to another or even within a state from one region to another.. Therefore, the focus of this study is to develop and customize accurate semi-empirical models to predict freezing and thawing depths in soils. The input data for these models must be readily available or can be obtained at a minimum cost. The other complex aspect of freeze-thaw cycles is the estimation of the frost heave of multi layered system due to ground freezing. During freezing, wet soils undergo heave due to the formation and growth of ice lenses. Frost heave is a function of many variables including soil type, its water holding capacity, and its thermal conductivity, air temperature, and frost depth. The heave could result in significant vertical and lateral stresses and movements which could lift 1 foundations or apply substantial additional stresses to exposed retaining structures. In addition, frost heave is typically followed by thaw consolidation and settlement. Therefore, spread footings located on soils subjected to freeze-thaw cycles would experience up and down vertical movements. For bridges, such movements may create unsafe driving conditions at the boundaries between the bridge and the adjacent pavement structure. On the other hand, as the spring begins the pavement starts to thaw from the top down and to a lesser extend from the bottom up (Marquis, 2008). The extend of frost and the time in which the pavement starts to thaw is a function of the material types, their thermal properties, water content, and climatic condition, such as temperature, wind speed, precipitation and solar radiation. At the beginning of the spring season , the pavement is in a critical condition where the upper and lower layers are thawed but the layer in between, which is still frozen, acts as an impermeable layer and trap water between the pavement surface and the undelaying materials causing saturation. As a results, the stiffness of the saturated layer and its bearing capacity decrease considerably leadingto substantial pavement deformations (Chapin et.al, 2012). Studies have shown that up to 90% of pavement damage occurs at this critical state (Tighe et.al, 2006). This phenomenon occurs particularly in low volume roads, these roads are often built by tight budgets and therefore they have minimal subbase and surface treatments. (Tighe et.al, 2006; Chapin et.al, 2012). Spring load restrictions (SLR) signs are usually placed along the roads as preservation strategies. The accuracy of the SLR implementation is critical in avoiding pavement damage. Even few days could lead to substantial damages. Further, the trucking industry should receive sufficient advance notice (at least 7-day notice) prior to posting the SLR in order to be prepared to follow the weight restriction. It is estimated that accurate posting and removing the SLR increase the life of low volume asphalt road’s by about 10 percent, which leads to a 2 potential saving of about $10,000,000 annually (Ovik et.al, 2000). Therefore an SLR policy which is simple, effective and accurate is essential. In this study, a simplified methodology was developed and verified to accurately estimate when to enforce SLR. 1.2 Study Objectives The objectives of this study are: 1. Review the advantages and shortcomings of some of the exsiting frost and thaw depth models. 2. Develop accurate and reliable models for predicting the frost depth during freezing period in Michigan. 3. Develop a model to predict heave and the resulting pressure under the pavement or behind existing retaining structures due to freezing of frost-susceptible soils in Michigan. 4. Develop accurate and reliable model for predicting the thaw depth under the pavement in spring season in Michigan. 5. Investigate changes in pavement bearing capacity in the cycles of freeze and thaw. 6. Develop a model to estimate when to post and remove SLR signs. 1.3 Research Plan To accomplish the objectives of the study, a research plan consisting of 4 tasks was drawn . The 4 tasks are summarized in the next few subsections. 1.3.1 Task 1 - Conduct Comprehensive Literature Review The literature review includes: 1. The state of the art of modeling of freeze and thaw in soils and their applicability to this study. 2. The state of art of modeling of frost heave and the resulting pressure behind retaining walls. 3 3. The state of the practice of SHAs for forecasting frost and thaw depth and the time at which to post and remove SLR . 1.3.2 Task 2 - Development of Heat Transfer Predictive Model After reviewing available models to predict the propagation of the freezing and thawing front in multilayered soils; the following steps will be taken: 1. The existing frost depth models that simulate the Michigan Department of Transportation (MODT) data the most will be further scrutinized and modified. 2. The thawing front will be modeled using a modified version of Nixon and McRoberts (1973) equation and other empirical models to fit the MDOT measured field data. 3. The models will be validated using Minnesota Department of Transportation (MNDOT) data. 1.3.3 Task 3 – Development of Coupled Heat and Mass Transfer Models for Prediction of Frost Heave and Frost Pressure The initiation and growth of ice lenses in a soil deposit in cold environment exert uplift pressure against the foundation and lateral pressures against a retaining structure and the soil behind it. Both pressures can be estimated using existing theory of coupled mass and heat transfer for estimating the rate of ice growth.The efforts in this task consist of the four steps listed below: 1. Estimate the freezing depth (Task 2). 2. Estimate the rate of flow of water to the frozen depth from a water supply (ground water table, surface water source, etc.) to calculate the rate of growth of ice lenses during the critical time period and at the most critical location of the site. The estimation of the rate of flow of water to the frozen depth could be based on several parameters including sub-zero 4 temperature, existence of frost susceptible soil, and the depth to the water table or the distance to the closest free supply of water. 3. Estimate the amount of heave pressure based on the frost heave model results. 4. Evaluate the accuracy of frost heave model using MDOT heave data. 1.3.4 Task 4 – Development of SLR Policy As ice in the pavement melts, the base and subgrade become saturated. Melting water is trapped in the upper subgrade, could not be drained through the frozen zone and consequently these layers lose strength. Since the frost-susceptible soils have relatively low hydraulic conductivity, it takes a long time for water to drain from the saturated layers. Therefore pavement weakness could continue for weeks after it is completely thawed. It is critical to implement and remove the load restrictions accurately. Therefore, the efforts in this task consist of the steps listed below: 1. The beginning of SLR period will be estimated using the surface temperature and thaw depth model. 2. The result of thaw depth model will be used to use to develop recommendations for removing the SLR signs. 1.4 Dissertation Layout This dissertation organized in 5 chapters and appendices. The contents of each chapter are detailed in the table of contents. The title of each chapter is listed below. Chapter 1 – Introduction & Research Plan Chapter 2 – Literature Review Chapter 3 – Data Mining Chapter 4 – Data Analyses & Discussion 5 Chapter 5 – Summary, Conclusions & Recommendations Appendices – Additional data, figures and drawings are presented in the appendices. 6 CHAPTER 2 LITERATURE REVIEW 2.1 Frost Depth One of the most important aspects of infrastructure design such as pavement, foundations, or utility line is frost depth prediction. Frost depth is a function of the material type, soil thermal properties, soil water content, and climatic conditions such as temperature, wind speed, precipitation, and solar radiation. In order to neutralize the effects of frost, foundations are usually built below the frost line. For pavements, most State Highway Agencies (SHAs) use nonfrost susceptible soils (granular materials). However, over time, fine aggregates migrate from the lower soil layers and soil becomes frost susceptible. In general, any soil might be considered frost susceptible when the percent fine (passing sieve number 200) exceeds about seven percent. Since silt has high water holding capacity and relatively low permeability, it is the most frost susceptible soil. Depending on the availability of the input data and the required accuracy, frost depth can be estimated by numerical, empirical, and/or mechanistic-empirical models. Numerical Models Different numerical techniques (finite element and finite difference) have been used for modeling complex transient heat flow in pavement layers. Hsieh et al (Hsieh et al., 1989) developed a three-dimensional finite difference computer program for predicting temperature profile in concrete pavements and rainfall infiltration into the layered system. The program inputs consist of typical meteorological year (TMY) data and typical physic al soil and concrete properties. They reported that their results were in a good agreement with the test results provided by the Florida DOT. 7 Dempsey (Dempsey 1985) developed a transient one-dimensional finite difference heat transfer model (Climatic-Materials-Structural (CMS)) to predict the temperature in asphalt layer. The required inputs for this model are thermal properties of materials, air temperature, solar radiation, and wind velocity. In this model, for temperature prediction the Fourier’s law was used as follows: 𝜕 2𝑇 𝜕𝑇 = 2 𝜕𝑧 𝛼𝜕𝑡 Equation 2.1 Where T = temperature (oC); z = depth (m); t = time (s); and α = thermal diffusivity (m2/s). The model has two boundaries, the surface temperature, and the temperature at the top of the base layer (McCartney et al., 2010). As can be seen in Figure 2.1, for estimating the surface temperature, the model considers heat convention and heat radiation in the energy balance equation as follows (ARA Inc., 2004) : 𝑄𝑖 − 𝑄𝑟 + 𝑄𝑎 − 𝑄𝑒 ± 𝑄𝑐 ± 𝑄ℎ ± 𝑄𝑔 = 0 𝑄𝑠 = 𝑄𝑖 − 𝑄𝑟 = 𝑎𝑠 𝑅 ∗ [𝐴 + 𝐵 𝑆𝑐 ] 100 Equation 2.2 Equation 2.3 9 0.28 𝑁𝑊 𝑄𝑎 = 𝜎𝑠𝑏 (𝑇𝑎𝑖𝑟 + 273.15) ∗ (0.77 − 0.074𝑝 ) (1 − ) 5 10 100 Equation 2.4 9 𝑁𝑊 𝑄𝑎 = 𝜎𝑠𝑏 𝜀((𝑇𝑠 + 273.15) ∗ )4 (1 − ) 5 100 Equation 2.5 𝑄𝑐 = ℎ (𝑇𝑎𝑖𝑟 − 𝑇𝑠 ) Equation 2. 6 ℎ𝑐 = 122.93 [0.00144 (𝑇𝑚 + 273.15)0.3 𝑈 0.7 − 0.00097(𝑇𝑠 − 𝑇𝑎𝑖𝑟 )0.3 ] Equation 2.7 8 Where Qi = incoming short wave radiation (W/m2); Qr = reflected short wave radiation (W/m2); Qa = incoming long wave radiation (W/m2); Qe = outgoing long wave radiation (W/m2); Qc = energy transferred to or from the body as a result of convection (W/m2); Qh = effect of transpiration, condensation, evaporation and sublimation (W/m2); Qg= energy absorbed by the ground (W/m2); Qs = net short wave radiation (W/m2); Ql = net long wave radiation (W/m2); as = surface short wave absorptivity of pavement surface; R* = extraterrestrial radiation incident on a horizontal surface at the outer atmosphere; A,B = constants that account for diffuse scattering and adsorption by the atmosphere; Sc = sunshine percentage; N = cloud base factor; W = average cloud cover during the day or night; Tair = air temperature (oC); σsb = Stefan-Boltzmann constant, 0.98 *10-8 (W/(m2.oC)); p = vapor pressure of the air (1 to 10mm Hg); ε = emissivity of the pavement; Ts = surface temperature (oC); hc = convection heat transfer coefficient; Tm = average of surface and air temperature (oC); and U = average daily wind speed (m/s). 9 Furthermore, this model was implemented by the Federal Highway Agency (FHWA) Integrated Climate Model (ICM) to investigate the environmental effect on the pavement. The ICM integrates the Infiltration and Drainage model (ID model), developed at Texas A&M University, the CMS model, the Frost Heave and Thaw Settlement model (CRREL model), developed at the United States Army Cold Region Research and Engineering Laboratory (CRREL) (Solaimanian and Bolzan, 1993, Johanneck, 2011). The ICM model considers that thermal properties of pavement layer do not change over time but for unbound layer the thermal properties vary due to the change in water and ice contents. Therefore, the pavement temperature predictions are coupled with the moisture estimations. As state before, the CMS is used to estimate the temperature within the asphalt layer. The boundary conditions for CRREL model are the ID and CMS outputs. Also, the required input variables are soil thermal conductivity, soil specific capacity, soil hydraulic conductivity and the soil-water characteristic curve (SWCC). The CRREL model is used to estimate the temperature within the subbase and subgrade layer. The model applies finite element solution to the governing water and temperature equations. The following equations are used for predicting the distribution of total hydraulic head and temperature, respectively (McCartney et al., 2010): 𝑑 𝑑𝐻 𝑑𝐻 (𝑘𝑢 )−𝑆( ) = 0 𝑑𝑧 𝑑𝑧 𝑑𝑡 Equation 2.8 𝑑 𝑑𝑇 𝑑 𝑑𝑇 (𝐶𝑤 𝑇) − 𝐶𝑠 ( ) = 0 (𝑘ℎ ) − 𝑑𝑧 𝑑𝑧 𝑑𝑧 𝑑𝑡 Equation 2.9 Where Ku = unsaturated permeability (m/s); H = total hydraulic head (m); S = slope of the soil water retention; 10 Kh = soil thermal conductivity (W/(m.oC)); Cw = water heat capacity (J/(Kg.oC)); Cs = soil heat capacity (J/(Kg.oC)); and All other parameters are the same as before. Figure 2.1 Heat transfer between pavement surface and air temperature (After Dempsey and Pur, 1990) 11 Solaimanian and Bolzan investigated the capability of the model in predicting the pavement temperature profiles accurately (Solaimanian and Bolzan, 1993). They performed a sensitivity analysis in order to evaluate the effect of different input parameters in the prediction results. The sensitivity analysis showed the following results: 1. While the air temperature significantly affects the pavement temperature predictions, the difference between air temperature and pavement temperature can be as low as 5.5 to 8 oC or as high as 22 to 28 oC depends on the solar radiation and percent sunshine. 2. In the same solar radiation and percent sunshine, air temperature and surface temperature has a linear relationship. 3. An increase in solar radiation from 395 to 526 (W/(m2)) leads to an increase of 4.5-5.5 oC in pavement temperature. 4. An increase in percent sunshine from 45% to 90% increases the pavement temperature by 4.5-5.5 oC in pavement temperature. 5. An increase in absorptivity from 0.7 to 0.8 or from 0.8 to 0.9 yields a 3 oC drop in pavement temperature at any depth. 6. An increase in emissivity from 0.7 to 0.8 or from 0.8 to 0.9 results in an increase of 3 oC in pavement temperature at any depth. 7. An increase in thermal conductivity from 1.7 to 3.4 and from 3.4 to 5.11 (W/(m.oC)) leads to a reduction of 2 oC and 1 oC in pavement surface temperature, respectively. However, the effect of thermal conductivity changes is found to be greater in larger depths. Solaimanian and Bolzan results showed that if the proper input variables were chosen the pavement temperature predictions were within ±1 oC of the measured surface temperature. 12 However, they recommended modification in user interface of the model in order to reduce the large number of required inputs variables (Solaimanian and Bolzan, 1993). Furthermore, in a collaboration study between the University of Illinois and Applied Research Associates (ARA), the ICM was modified and its moisture prediction capability was improved. The modified model was called Enhanced Integrated Climate Model (EICM) (Zapata and Houston, 2008). Since then, the accuracy of the results of EICM was investigated by different researchers. Liang (Liang, 2006) used EICM to estimate the temperature and moisture profile in six pavement section in the State of Ohio. The results showed that the predicted temperature profiles did not coincide with the measured data but they could be considered within an acceptable range. They also found that there is a good correlation between frost depth predictions and measured ones in sections with unbound base materials but not in sections with bounded base material. Heydinger (Heydinger, 2003) evaluated the EICM temperature prediction in two sites in the State of Ohio for the year 2000. He used the default input data in the model and his results showed that the EICM consistently over predicted pavement temperature. Ahmed et al (Ahmed et al, 2005) used site specific input values and compared the EICM temperature and moisture predictions with the measured values in different site in New Jersey. Their results showed that while the predicted temperatures follow the same trend as the measured ones, their difference can be significant. However, other researchers’ results indicated that after local calibration, EICM temperature predictions are relatively accurate (Khazanovich, 2013; Chung, and Shin, 2015). Yavuzturk et al (Yavuzturk et al., 2005) proposed a transient two-dimensional finite difference model to assess the thermal behavior and temperature distribution in asphalt pavement. TMY weather data were used and sensitivity analyses were conducted to determine 13 the influence of different thermal properties of the materials on the predicted asphalt temperature. They reported that the temperature predictions were most effected by variation of the absorptivity, volumetric heat capacity, emissivity, and thermal conductivity of the materials. Chapin et al (Chapin et al., 2012) utilized finite element program TEMP/W (GEO-SLOPE 2007) to simulate freezing and thawing front in the pavement. They applied the program to two sites in northern Ontario with considerably different pavement structures. First, a steady state analysis was conducted to establish the initial conditions within the model and second, a transient analysis was conducted. By using adiabatic1 conditions on the lateral boundaries they induced one-dimensional heat flow. They reported that the predicted frost front was several days behind the measured frost front. 2.2.2.1 UNSAT-H Modeling UNSAT-H is a one dimensional, finite difference computer program developed at the Pacific Northwest Laboratory (Fayer and Jones, 2000). UNSAT-H can simulate the water and heat balance in a layered cross section simultaneously. The input properties for the models are listed below: 1. Hydraulic Properties - To solve the water balance equations, relationships for both water content and hydraulic conductivity as a function of suction head are required. To describe soil water retention from measured data the van Genuchten function has been used: 1 An ADIABATIC process is the changing temperature of air due to its movement. Rising air will cool adiabatically, whereas sinking air warms adiabatically. The DIABATIC process, on the other hand, is any change in air temperature not associated with adiabatic vertical displacement of air. The prime source of heating in the DIABATIC process is the sun, while the main cause of cooling is evaporation and the emission of long wave energy from the ground surface. 14 𝜃 = 𝜃𝑟 + (𝜃𝑠 − 𝜃𝑟 )[1 + (𝛽ℎ)𝑛 ] −𝑚 Equation 2.10 Where 𝜃𝑟 = residual water content; 𝜃𝑠 = saturated water content; h = suction; and 𝛽, n, m = fitting parameters 2. Thermal Properties - UNSAT-H model use Cass et al. equation to express thermal conductivity (𝑘ℎ ) as a function of water content (Stormont and Zhou, 2001): 𝜃 𝜃 −[𝐶 ] 𝑘ℎ = 𝐴 + 𝐵 − (𝐴 − 𝐷) exp 𝜃𝑠 𝜃𝑠 𝐸 Equation 2.11 Where 𝑘ℎ = thermal conductivity (W/(m.oC); 𝜃 = the water content corresponding to the measured 𝑘ℎ ; A, B, C ,D, E = the fitting parameters; and All other parameters are the same as before. Mechanistic Empirical Models Neumann proposed the first solution to the heat transfer phase-change problem in his lectures in the 1860’s; he then published his work in 1912 (Jiji, 2009). In his solution, onedimensional heat transfer in a semi-infinite region was assumed. The above freezing initial surface temperature (Ti) drops to T0 (a temperature below the freezing point) and freezing starts to propagate through the liquid phase as shown in Figure 2.2 (Jiji, 2009). The governing heat conduction equations for solid and liquid phases are stated in Equation 2.12 and Equation 2.13, respectively. 𝜕 2 𝑇𝑓 1 𝜕𝑇𝑓 = 2 𝜕𝑥 𝛼𝑓 𝜕𝑡 15 0<𝑥<𝑃 Equation 2.12 𝜕 2 𝑇𝑢 1 𝜕𝑇𝑢 = 𝜕𝑥 2 𝛼𝑢 𝜕𝑡 𝑥>𝑃 Equation 2.13 Where the subscripts u and f refer to unfrozen and frozen, respectively; t = time since the freezing starts (s); 𝑥𝑖 = frost depth (m); T = temperature (oC); and 𝛼 = thermal diffusivity (m2/s) calculated using Equation 2.14. 𝛼= 𝑘 𝜌𝑐𝑝 Equation 2.14 Where 𝑘 = thermal conductivity of the soil (W/(m.oC)); 𝑐𝑝 = specific heat at constant pressure (J/(Kg.oC)); 𝜌 = density (Kg/m3). Moving Interface Frozen zone Unfrozen zone xTf (x,t) P Figure 2.2 A schematic representation of two-phase heat conduction The interface energy equation is stated in Equation 2.15: 𝜕 2 𝑇𝑓 (𝑥𝑖 , 𝑡) 𝜕 2 𝑇𝑢 (𝑥𝑖 , 𝑡) 𝑑𝑥𝑖 𝑘𝑓 − 𝑘𝑢 = 𝜌𝑓 𝑙 2 2 𝜕𝑥 𝜕𝑥 𝑑𝑡 The boundary conditions are 𝑇𝑓 (0, 𝑡) = 𝑇0 16 Equation 2.15 𝑇𝑓 (𝑥𝑖 , 𝑡) = 𝑇𝑚 𝑇𝑢 (𝑥𝑖 , 𝑡) = 𝑇𝑚 𝑇𝑢 (∞, 𝑡) = 𝑇𝑖 And the initial conditions are 𝑇𝑢 (𝑥, 𝑡) = 𝑇𝑖 𝑥𝑖 (0) = 0 Where the subscripts u and f refer to unfrozen and frozen, respectively l = latent heat of fusion (J/Kg); 𝑇𝑚 = bulk freezing temperature (oC); and All other parameters are the same as before. The frost depth can be estimated using Equation 2.16: 𝑃 = 𝜇 √4𝛼𝑓 𝑡 Equation 2.16 Where P = frost depth (m); 𝜇 = constant obtained from; and All other parameters are the same as before. The parameters 𝜇 can be calculated using Equation 2.17: 𝛼𝑓 𝑘𝑢 𝑇𝑖 − 𝑇𝑚 exp(−𝜇2 ) −√ 𝑒𝑟𝑓𝜇 𝛼𝑢 𝑘𝑓 𝑇𝑚 − 𝑇0 𝜇 2 𝛼𝑓 exp (− 𝛼 ) 𝑢 𝛼𝑓 1 − erf (√𝛼 𝜇) 𝑢 = √𝜋𝜇𝑙 𝑐𝑝𝑓 (𝑇𝑚 − 𝑇0 ) Where the subscripts u and f refer to unfrozen and frozen, respectively; erf = Gauss error function; Ti = initial surface temperature (oC); T0 = surface temperature at t≠0 (oC); and 17 Equation 2.17 All other parameters are the same as before. Further, Stefan solved Neumann’s equation for a special case of no heat transfer in liquid layer in 1891, (Jiji, 2009) as follow: 𝑃=√ 2𝑘𝑓 (𝑇𝑚 − 𝑇0 ) 𝑡 𝜌𝑙 Equation 2.18 Where all parameters are the same as before. Stefan assumed that the applied constant surface temperature (𝑇𝑚 − 𝑇0 ) multiplied by the time (t) is equivalent to the freezing index (FI) at that time. He further introduced a dimensionless multiplication parameter (n) to converts air temperature to surface temperature. then , Equation 2.19 became: 172.8 𝑘𝑓 ∗ 𝑛 ∗ 𝐹𝐼 𝑃=√ 𝐿 Equation 2.19 𝐿 = 334𝑤𝛾𝑑 Equation 2.20 Where P = depth of freeze or thaw (m); 𝑘𝑓 = thermal conductivity of soil (W/(m.oC)); n = dimensionless parameter which converts air index to surface index; FI = freezing index (oC-day); note that the freezing index in Stefan equation is similar to the cumulative degree-day at time t, it is not the conventionally defined freezing index for a winter season. L = volumetric latent heat of fusion (KJ/m3); w = water content and; and 𝛾𝑑 = dry density (Kg/m3). 18 Since Stefan’s equation does not consider the volumetric heat capacity of the soil and water the accuracy of the results are debatable. Consequently, several studies have been conducted to improve the prediction of frost depth, including the modified Berggren’s equation (Aldrich et al., 1953). Berggren’s Equation is very much similar to the early work of Neumann. Therefore, it is not explained here. Aldrich et al. applied a correction factor to Berggren’s equation, which is a function of two dimensionless parameters, the thermal ratio (α), and the fusion parameter (μ) (see Figure 2.3). In this figure V0 is the initial temperature differential (mean annual temperature -0 0C), Vs is the average temperature differential (nFI/t), C is the average volumetric heat capacity, and L is the volumetric heat of fusion. These parameters take the effect of temperature changes in the soil mass into account and depend on the freezing index, the annual average temperature in the site and the thermal properties of the soil (USACE, 1988). The modified Berggren’s Equation can be written as follows: 𝑃 = 𝜆√ 172.8 𝑘 ∗ 𝑛 ∗ 𝐹𝐼 𝐿 Equation 2.21 Where 𝜆 = correction factor; and All other parameters are the same as before. A multilayer solution to the modified Berggren’s equation can be applied to nonhomogeneous soils by calculating the required cumulative freezing degree-day (CFDD) for frost to penetrate each layer. The maximum summation of the CFDDs must be equal to or less than the regional and seasonal freezing index. The frost depth can be estimated as the sum of the thicknesses of all the frozen layers (USACE, 1988). The CFDD required to penetrate the nth layer is defined as: 19 𝐶𝐹𝐷𝐷𝑛 = 𝑛−1 𝐿𝑛 𝑑𝑛 𝑅𝑛 [(∑ 𝑅𝑛 ) + ] 2 86.4𝜆𝑛 2 1 Equation 2.22 Where Ln = volumetric latent heat of fusion of the nth layer (KJ/m3); Rn= thermal diffusivity of the nth layer= dn/kn ((m2.oC)/W); dn = depth of the nth layer (m); 𝜆𝑛 = correction factor of the nth layer; kn = thermal conductivity of the nth layer (W/(m.oC)); and CFDDn = cumulative freezing degree day required for frost to penetrate the nth layer (°Cdays). Figure 2.3 Fusion parameter (μ) versus correction factor (λ) 20 The Pavement-Transportation Computer Assisted Structural Engineering (PCASE) software provided a more accurate numerical solution of the Modified Berggren’s equation, (Bianchini et al. 2012). Berg (Berg, 1996) applied Modified Berggren’s equation to 40 sites in the state of Minnesota for 3 years to assess the accuracy of the results. He reported that predicted frost depths were within ±15 percent of the measured frost depth. He also conducted different sensitivity analysis to assess the dependence of the predicted frost depths to the n-factor (defined on page 2-5), water content, dry density, thermal conductivity, and each layer thickness. Berg concluded that small variation in thickness, water content and dry density of each layer would have a small effect on the predicted frost depths. On the other hand he found that increases in the n-factor values would result in deeper frost depths prediction. Whereas increasing the measured thermal conductivity by 25 percent would lead to better frost depths prediction. Stated differently, Berg found that the modified Berggren’s equation produced more accurate estimates of the frost depth when the measured thermal conductivity was artificially increased by 25 percent. Empirical Models Chisholm and Phang used the data from different stations throughout Ontario and developed an empirical equation to correlate the calculated cumulative freezing degree day (CFDD) and the measured frost depths (Chisholm and Phang, 1983). 𝑃 = 0.0578 √𝐶𝐹𝐷𝐷 − 0.328 Where P = depth of freeze or thaw (m); and All other parameters are the same as before. 21 Equation 2.23 Many State Highway Agencies (SHAs) used similar approach to generate their own equations or simply calibrated Equation 2.23 using local frost depth data and CFDD. Dore (Tighe et al, 2007) conducted a research to develop an empirical model for frost depth in Quebec, Canada. First, he developed Equation 2.24 to estimate pavement surface temperatures (PST) based on the measured air temperatures. Second, he calculated the cumulative freezing degree day (CFDD) based on the estimated pavement surface temperature (PST of Equation 2.24) and estimated the frost depth using Equation 2.25. Third, he correlated the estimated frost depths from Equation 2.25 to the measured frost depth and obtained statistical Equation 2.26. 𝑃𝑆𝑇 = 𝑇𝑀𝐸𝐴𝑁 + [0.178(𝑇𝑀𝐴𝑋 − 𝑇𝑀𝐼𝑁 )] + 1.628 Equation 2.24 𝑃 = 𝐶 √𝐶𝐹𝐷𝐷 Equation 2.25 2 0.5 𝑃𝑐𝑜𝑟𝑟 (√𝐶𝐹𝐷𝐷 − 𝑋𝑀𝐸𝐴𝑁 ) 1 = 𝑃 + [𝐶𝐼(𝑆𝑒 ) (1 + )+( ) ∑(𝑋𝑖 − 𝑋𝑀𝐸𝐴𝑁 )2 398 ] Equation 2.26 Where; TMEAN = (TMAX+TMIN)/2; TMAX = maximum daily air temperature (oC); TMIN = minimum daily air temperature (oC) and; and PST = estimated pavement surface temperature (oC). P = frost depth (cm); C = regression constant; and CFDD = cumulative freezing degree days based on the estimated pavement surface temperature (PST) (oC-day). Pcorr = corrected frost depth; 22 CI = confidence interval for a population mean, a function of significance level, alpha = 0.4, one standard deviation and a sample size of one; Se = sum of squared errors; Xi = measured frost depth (cm); and XMEAN = Average measured frost depth (cm). Tighe et.al (Tighe et.al, 2007) used data from one study site along Highway 569 in Northern Ontario and calibrated the Chisholm and Phang model. Furthermore, they used CFDD and cumulative thawing degree day (CTDD) and developed a modified model for estimating the frost depths as follow: For 0 ≤ 𝑖 ≤ 𝑖0 For i ≥ 𝑖0 𝑃𝑖 = 𝑎 + 𝑏√𝐶𝐹𝐷𝐷𝑖 + 𝑐√𝐶𝑇𝐷𝐷𝑖 Equation 2.27 𝑃𝑖 = 𝑑 + 𝑒√𝐶𝐹𝐷𝐷𝑖 + 𝑓√𝐶𝑇𝐷𝐷𝑖 Equation 2.28 Where i = number of days after the day indexed as day i=0 (i= 0 day on which air temperature first falls below 0 oC); i0 = day after which the CTDD consistently increases; Pi = depth of frost on day i; CFDDi = cumulative freezing degree day on day i (oC-days); CTDDi = cumulative thawing degree day on day i (oC-days); and a,b,c,d,e,f = calibration coefficients. Moreover, they used Road Weather Information Systems (RWIS) in three sites close to the study site to estimate the frost depths and compare them with the study site data. They estimated the calibration coefficients and calculated the frost depth. Although, the coefficient of determination was 91%, the reliability of the model is questionable since only one year of data was used. 23 2.2 Frost Heave In seasonally frozen regions, soil freezing causes frost heave, which may cause extensive damage to various civil engineering structures, such as pavements and utility lines (Liu et.al, 2013). Frost heave refers to the uplifting of ground surface caused by freezing of water within the layers of soil. Taber (Taber, 1930) was the first one that demonstrated experimentally the features of frost heave. Before Taber frost heaving was explained based on experiments with closed systems. Taber showed that under normal conditions, the freezing occurs in an open system. Therefore in the freezing process water migrates through the soil voids below the freezing zone, causes excessive heaving by creating segregated ice layers. Tendency of a soil to heave under the freezing conditions is known to be influenced by parameters such as soil type, freezing rate, availability of water and the applied load or overburden pressure Frost Heave Mitigation The effects of frost heave on various structures vary from one structure to the next. Typically, structural foundations are constructed below the expected frost depths and hence, they are not affected by frost heave. Frost susceptible soils or free standing water behind bridge abutments and/or behind exposed retaining structures (such as retaining structures along depressed highways), are subjected to frost and frost heave causing active pressure against the structures. Basement retaining walls are rarely affected by frost due to heat loss from the basement interior that keeps the soil in the vicinity of the wall in relatively warm conditions. Pavement structures are frost heave susceptible especially if the roadbed soils are not protected from frost action or if the granular base and subbase are subjected to saturation due to lack of proper drainage. Given the potential damage due to frost heave, different techniques have been 24 proposed to mitigate frost heave damage especially in the pavement. The most common techniques are: 1. Cutting off the Water Source - The source of water can be cut off in many different ways. One common technique is to install a barrier between the water source and the frost zone (Edgar, 2014). The barrier reduces the capillary action and consequently reduces frost heave. A blanket or a layer of gravel and crushed stone under the pavement or wrapping the roadbed soil by a geo-membrane layer could be effective in decreasing access to water (Wallace, 1987; Edgar, 2014). Another technique is to remove water using a proper drainage system. In pavements, drain tile, edge drain, and/or open side ditches can be built to remove the water. In retaining wall, weep holes can be installed at the foot of the wall which is exposed to frost (Wallace, 1987). 2. Removing Frost Susceptible Soil - As stated in the previous section, some soils are more frost susceptible than others. Such soils can be replaced by non-susceptible soils if the cost is not prohibitive. In a typical scenario, the various frost heave mitigation options are assessed against their costs. The most cost effective option is typically chosen. 3. Reducing Freezing Depth – Although different approaches can be used to prevent frost penetration, two of these approaches are insulation and chemical additives to lower the water freezing temperature. Since insulation is the most common method, it is detailed further below. 4. Insulation Method – This method could be used in many different structures to decrease heat loss from the soil to the atmosphere. In pavements, an insulation layer is typically placed above the roadbed soils to protect the soils from freezing. Rigid polystyrene foams (RPF) are commonly used for frost protection under different building foundation and infrastructures. 25 Two types of polystyrene have been used; expanded polystyrene (EPS) and extruded polystyrene (XPS). The insulation materials are usually known by their thermal resistivity (R-value). R-value is an indication of material resistance to heat flow. It has inverse relationship with thermal conductivity of the material (Edgar, 2014). Table 2.1 shows the R-value of different RPF according to ASTM C578. It should be noted that the nominal R-value varies depending on moisture exposures condition. Moisture condition could vary from one site to another and it depends on the drainage system and on the direction along which the insulation is installed (vertical or horizontal). Therefore in the design process, the effective R-values are calculated or estimated and used. Table 2.1 Thermal resistance values (R-values) at different mean temperature Classification XI I VII II IX XIV Minimum Density 11.2 14.4 18.4 21.6 28.8 38.4 (kg/m3) Mean Temperature -3.9 ± 1 oC 4.4 ± 1 oC 43.3 ± 1 oC XII X IV VI VII V 19.2 20.8 23.2 28.8 35.2 48.1 Thermal Resistance of 2.54 Centimeters Thickness Minimum ((m2.oC)/W) 1.97 2.39 2.51 2.62 2.74 1.88 2.28 2.39 2.51 2.62 1.65 1.85 1.97 2.08 2.19 2.74 2.85 2.19 2.96 3.08 2.45 3.19 3.08 2.65 3.19 3.08 2.65 3.19 3.08 2.65 3.19 3.08 2.65 3.19 3.08 2.65 Another important property of the PRF is the minimum thickness. Non-uniform distribution of moisture in RPF leads to edge effects and as the insulation thickness decreases it impacts the thermal performance of the RPF. It should be noted that the effect of the thickness varies depending on the insulation type and moisture conditions (Crandell, 2010). Table 2.2 shows the Design values for frost protected shallow foundation (FPSF) RPF based on ASCE 3201. 26 Various researches investigated the effect of temperature and moisture conditions on the RPF properties. Ojanen and Kokko (1997) used data from different highway projects to evaluate the EPS performance. They found that the thermal conductivity measured at -5 oC is the most relevant to the highway conditions. Their data showed that with proper drainage, the long term moisture contents in EPS under highways are in the range of 0.5 to 2.5 %. Sandberg (1986) did research on RPF performance under highways. He found that moisture content distribution is highly non-uniform in XPS which reduces the influence of moisture content on R-value in comparison to EPS (Crandell, 2010). Nevertheless, it is apparent that XPS performs consistently under different conditions. But EPS performance could vary based on the moisture content, density and manufacturing process (Crandell, 2010). Table 2.2 Design values for FPSF insulation materials based on ACSE 32-01 Insulation Type Minimum Density (kg/m3) Effective Resistivity (R/cm.) Vertical Horizontal Nominal Resistivity (R/cm.) Allowable Bearing Capacity (kg/m2) Minimum Insulation Thickness (cm.) Vertical Horizontal ESP II IX 21.6 28.8 0.7 0.8 0.6 0.6 0.9 0.9 N/A 5859 5.1 3.8 7.6 5.1 1.1 1.1 1.1 1.1 1.1 N/A 5859 9374 14061 23436 3.8 2.5 2.5 2.5 2.5 5.1 3.8 2.5 2.5 2.5 XPS X IV VI VII V 21.6 25.6 28.8 35.2 48.1 1.0 1.0 1.0 1.0 1.0 0.9 0.9 0.9 0.9 0.9 The magnitude and the rate of frost heave can be predicted in terms of certain characteristics of the freezing system and some boundary conditions by use of a practical theory explaining the frost heave of a specific soil (Konrad and Morgenstren, 1980). In general the theories toward this matter can be classified into two categories, capillary theory and frozenfringe theory. 27 2.2.2.1 Capillary Theory Capillary theory, also known as primary frost heave theory, is characterized by a frozen and an unfrozen zone within the soil strata. Consider pure water to be at equilibrium with ice, when a differential amount of water freezes at constant temperature and pressure: For two phases of ice and water; 𝑑𝐺 = 𝑉𝑑𝑃 − 𝑆𝑑𝑇 = 0 Equation 2.29 𝑑𝐺𝑖 = 𝑑𝐺𝑤 Equation 2.30 𝑉𝑖 𝑑𝑃 − 𝑆𝑖 𝑑𝑇 = 𝑉𝑤 𝑑𝑃 − 𝑆𝑤 𝑑𝑇 Equation 2.31 By rearrangement the equation becomes 𝑑𝑃 𝑆𝑖 − 𝑆𝑤 ∆𝑆𝑤𝑖 = = 𝑑𝑇 𝑉𝑖 − 𝑉𝑤 ∆𝑉𝑤𝑖 Equation 2.32 Where the subscripts “i” and “w” stand for ice and water, respectively (Takagi, 1978). G = Gibbs free energy (J); S = entropy (J/oC); V = volume (m3); P = pressure (Pa); and T= temperature (oC). The entropy change ∆𝑆𝑤𝑖 and the volume change ∆𝑉𝑤𝑖 are the changes, which occur when a unit amount of water is transferred from phase w to phase i at the equilibrium temperature and pressure. Clapeyron substituted the latent heat of phase transition as ∆𝐻𝑤𝑖 = 𝑇𝑚 ∆𝑆𝑤𝑖 in Equation 2.32 and obtained Equation 2.33 (Smith et al. 2001): 𝑑𝑃 ∆𝐻𝑤𝑖 = 𝑑𝑇 (𝑇𝑚 + 273)∆𝑉𝑤𝑖 28 Equation 2.33 Where ∆𝐻𝑤𝑖 = the enthalpy change when a unit amount of water is transferred from water to ice; Tm = the bulk freezing temperature, (oC); and All other parameters are the same as before. Equation 2.33 can be rewrite as (Peppin and Style, 2012): 𝑃𝑖 − 𝑃𝑤 = 𝜌𝑤 𝐿 (𝑇 − 𝑇) 𝑇𝑚 𝑚 Equation 2.34 Where T = the thermodynamic equilibrium temperature of the system (oC); L = latent heat of fusion (J/kg); 𝑃𝑖 = ice pressure (Pa); 𝑃𝑤 = water pressure (Pa); and 𝜌𝑤 = density of water (kg/m3). The Clapeyron equation explains thermodynamically why lowering the temperature below the freezing temperature causes water to move (be sucked) toward the ice. Black (Black, 1995) solved Clapeyron equation for different scenarios. 1. If the pressure difference in ice and water are the same, by increasing the confining pressure of 1 MPa the melting temperature decreases by 0.074 oC. 2. If the change in confining pressure in water is 1.09 times greater than the change in ice pressure then the melting temperature remains constant. 3. If water pressure is constant by increasing the ice confining pressure of 1 MPa the melting temperature decreases by 0.893 oC. 4. If ice pressure is constant by decreasing the water confining pressure of 1 MPa the melting temperature decreases by 0.810 oC. Everret (Everret, 1960) constructed a simple model for explaining the capillary theory. He considered two cylinders each closed up by a piston and joined by a capillary tube as shown in 29 Figure 2.4. By lowering the temperature, water starts to freeze in the upper cylinder. When the upper cylinder is completely filled up with ice; further decreases in temperature result in water flow from the lower cylinder to the upper one. According to Laplace equation if the radius of the capillary tube is r (m), ice can only penetrate to the capillary tube when Equation 2.35 is satisfied. 𝑃𝑖 − 𝑃𝑤 = 2𝜎𝑖𝑤 𝑟 Equation 2.35 Where 𝜎𝑖𝑤 = ice-water surface energy (J/m2); All other parameters are the same as before. Ice Ice r=Radius of capillary tube Ice Pi Capillary tube r Pw p Figure 2.4 (a) water Figure 2.4 (b) Figure 2.4 a. Equilibrium between ice and water, b. simple ice-water model Since the capillary tube represents the soil pores, it implies that segregated ice forms when 𝑃𝑖 < 𝑃𝑤 + And pore ice forms when 30 2𝜎𝑖𝑤 𝑟 Equation 2.36 𝑃𝑖 > 𝑃𝑤 + 2𝜎𝑖𝑤 𝑟 Equation 2.37 This implies that the growth of ice lenses will stop as the ice invades the soil at the maximum heaving pressure given in Equation 2.38 (Loch and Miller, 1975): 𝑝𝑚𝑎𝑥 = 𝑃𝑤 + 2𝜎𝑖𝑤 𝑟 Equation 2.38 The temperature Tl at which ice invades the pores can be found by combining Equation 2.34 and Equation 2.35 into Equation 2.39 (Peppin and Style, 2012): 𝑇𝑙 = 𝑇𝑚 (1 − 2𝜎𝑖𝑤 ) 𝑟𝜌𝑤 𝐿 Equation 2.39 Where 𝑇𝑙 = The temperature at which ice invades the pores (oC), and All other parameters are the same as before. The capillary theory has various limitations including: 1. Predictions of the maximum frost-heave pressure works well with idealized soils composed of particles with one size. But, in soils with different particle sizes the heaving pressures are considerably larger (Peppin and Style, 2013). 2. Capillary theory can be used to predict the flow rate towards the ice lenses in the frozen region. By assuming that the porous medium is incompressible, Darcy’s law can be used to determine the flow rate of water towards the lenses by Equation 2.40. But the equation tends to over predict the measured values of flow rate (Peppin and Style, 2012) 𝑉= 𝑘 𝑃𝑅 − 𝑃𝑓 𝜇 𝑍ℎ Where V =flow rate (m/s); k = permeability of the soil (m/s); μ = dynamic viscosity of water (Pa.s); 31 Equation 2.40 Zh = distance between the ice lens and the water reservoir (m); PR = ground water pressure (Pa); and Pf = pressure of the water directly below the warmest lens (Pa). 1. No mechanism for initiation of new lenses has been explained by this method. 2.2.2.1 Frozen Fringe Theory Since the capillary theory has limitations, some researchers explained the propagation of frost heave phenomenon by another theory, Frozen Fringe theory. Frozen Fringe theory, also termed as secondary frost heave, is characterized by three zones: a frozen zone, a partially frozen zone and unfrozen zone. According to this theory, frost heave can continue to occur at ice-lens temperatures above 𝑇𝑙 (the temperature at the bottom of the frozen zone) when a frozen fringe is shaped by formation of ice in the soil pores, see Figure 2.5. Figure 2.5 Schematic diagrams for the frost heave process (Peppin and Style, 2012) Stated differently, if the rate of extracting heat is too large or the soil column is too tall, or too impervious to prevent ice entry the frozen fringe is created beneath ice lenses at the top. 32 Therefore ice pressure could rise above the maximum (Pmax).This process is called secondary heaving (Loch and Miller, 1972). At the interface of ice lens and soil particles, there are repulsive intermolecular forces (surface tension). These forces act like a disjointing pressure that separate ice and soil particles and initiating a microscopically thin layer of water between the ice lenses and the soil particles below the freezing temperature, Tm (Dash et al. 2006), see Figure 2.6. Because of the repulsive forces between the ice lenses and the soil particles, the pressure in the thin water film is reduced causing suction and upward water movement toward the growing ice lenses (Peppin and Style, 2012). Figure 2.6 Schematic diagram of a freezing soil with frozen fringe (After Peppin and Style, 2012) Secondary frost heave can be affected by the suction pressure. The specific characteristics of the soil determine the practical relation between the suction and the unfrozen water content. As the ice-water interface curvature is increasing, the unfrozen water content decreases which consequently yields an increase in suction. 33 According to the Clapeyron formula, an increase in load results into a decrease in the amount of unfrozen water, which consequently increases the suction. However, increase in load makes the onset of new ice lens formation more difficult. Therefore, higher suction is required to separate the soil grains. These dual effects of increase in load make the ice lenses initiation in clays more easily. The secondary frost heave occurs mostly under any load condition in clays, while rarely takes place in sands (Fowler and Krantz, 1994). When the freezing front penetrates into the soil, it absorbs the moisture in the soil, which stands for a process of both heat and mass transfer (Harlen1973). The complexity in the frost heave theory arose from this coupled effect of heat and mass transfer. The first model, which considers heat and mass flow in the soil, was proposed by Harlan (Harlan, 1973). He proposed that the generalized one-dimensional mass flow for steady or unsteady flow in a saturated or partially saturated soil media can be modeled by Equation 2.41 and the onedimensional transient heat transfer can be modeled by Equation 2.42 𝜕 𝜕𝐻 𝜕(𝜌𝑤 𝜃𝑙 ) [𝜌𝑤 𝐾(𝑥, 𝑇, 𝜓) ] = + ∆𝑀 𝜕𝑥 𝜕𝑥 𝜕𝑡 Equation 2.41 𝜕 𝜕𝑇 𝜕(𝜗𝑥 𝑇) 𝜕(𝐶𝑇) [𝑘(𝑥, 𝑇, 𝑡) − 𝑐𝑙 𝜌𝑤 = 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑡 Equation 2.42 Where t = time (min); ρw = density of water fraction (gr/cm3); θl = volumetric water content = volume of water/ total volume (cm3/cm3) ; K = hydraulic conductivity (cm/ min); T = temperature (Co); H = total head (cm); 𝜓 = capillary pressure head (cm); 34 ∆𝑀 = change in mass of ice per unit volume, unit time (gr/(cm3. min)); cl = bulk specific heat of water (cal /gm /Co); 𝜗𝑥 = water flow velocity in x direction (cm/ min); C = 'apparent' volumetric specific heat (cal /cm2/Co); and All other parameters are the same as before. Gilpin (Gilpin 1979, 1980a, 1980b) studied water flow towards the ice layer and proposed a physical model for prediction of ice lensing and heave rate, he suggested that frost heave is a function of basic soil properties and boundary conditions. He assumed that the free energy of water in the pores is lowered by the surface effect of the solid. Figure 2.7 shows the pressures in the water near the solid soil surface in the case of the existence of tension between the water meniscus and the soil. The effect of the tensile surface force on free energy could be described as follow: 𝐺𝑤 = 𝐺𝑤𝑂 + 𝜈𝑤 𝑃𝑤 − 𝑆𝑤 𝑇𝑤 − 𝑔(𝑦) Equation 2.43 Where the subscripts w stand for water; g(y) = a dummy variable expressing the effect of the particle surfaces on the free energy of water (KJ/mol), g(y) is estimated using Equation 2.44 by setting y equal to h (the distance between the soil particle surface and the ice lenses); 𝐺𝑤 = the free energy of water near the surface (KJ/mol); 𝐺𝑤𝑂 = the free energy at bulk conditions T0 and P0 (KJ/mol) ; νw = specific volumes of water (m3/Kg); and All other parameters are the same as before. 35 In the case of thermodynamic equilibrium, Gw should be constant in the water layer and also the temperature could be considered constant because the layer is thin. Therefore the effect of surface can be obtained as: 𝑃𝑤𝑦 = 1 𝑔(𝑦) 𝜈𝑤 Equation 2.44 Where 𝑃𝑤𝑦 = pressure at the distance y from the surface; and All other parameters are the same as before. Figure 2.7 A schematic representations of equilibrium conditions for ice and water near a substrate (Gilpin, 1980) On the other hand, the free energy in the ice (Gi) can be obtained as: 𝐺𝑖 = 𝐺𝑖𝑂 + 𝜈𝑖 𝑃𝑖 − 𝑆𝑖 𝑇𝑖 Where the subscripts i stand for ice 36 Equation 2.45 𝐺𝑖 = the free energy of ice (KJ/mol); 𝐺𝑖𝑜 = the free energy at bulk conditions T0 and P0 (KJ/mol); νi = specific volumes of ice (m3/Kg); and All other parameters are the same as before. The temperature cannot be different across the phase boundary, but the pressure difference can be gained as: ̅ 𝑃𝑖 − 𝑃𝑤ℎ = 𝜎𝑖𝑤 𝐾 Equation 2.46 ̅ = the mean curvature of the interface, Where 𝐾 𝑃𝑤ℎ = pressure at the distance h from the surface, and All other parameters are the same as before Equating Equation 2.43 and Equation 2.45 and using Equation 2.46, g(h) can be calculated as ̅− 𝑔(ℎ) = −Δ𝜈𝑃𝑤ℎ − 𝜈𝑖 𝜎𝑖𝑤 𝐾 𝐿𝑇 𝑇𝑚 Equation 2.47 Using Equation 2.46 and Equation 2.47, the pressure gradient can be obtained as 𝑑𝑃𝑤𝑦 𝜈𝑖 𝑑 𝐿𝑇 𝜈𝑖 𝑑 𝐿𝑇 ̅+ = [𝑃𝑤ℎ + 𝜎𝑖𝑤 𝐾 ]= [𝑃𝑖 + ] 𝑑𝑥 𝜈𝑤 𝑑𝑥 𝜈𝑖 𝑇𝑚 𝜈𝑤 𝑑𝑥 𝜈𝑖 𝑇𝑚 Equation 2.48 Therefore the flow rate through water layer (q) can be calculated as 𝑞 = −𝑘 𝜈𝑖 𝑑 𝐿𝑇 [𝑃𝑖 + ] 𝜈𝑤 𝑑𝑥 𝜈𝑖 𝑇𝑚 Equation 2.49 The above equation can be modified to obtain the water velocity in the frozen fringe as 𝑉𝑓𝑓 = 𝐾𝑓 𝜈𝑖 𝑑 𝐿𝑇 [𝑃𝑖 + ] 𝑔 𝑑𝑥 𝜈𝑖 𝑇𝑚 Where 𝐾𝑓 = the permeability in the frozen fringe (m/s), 𝑉𝑓𝑓 = the velocity of water flow in the frozen fringe (m/s); and 37 Equation 2.50 All other parameters are the same as before. Figure 2.8 illustrates schematically the frost heave simulation. A linear temperature profile is assumed in each layer and the heat balance equation can be written as: Figure 2.8 the frost heave simulation model (Gilpin, 1980) − 𝑘𝑓 (𝑇𝑇𝑂𝑃 − 𝑇𝑙 ) 𝑘𝑓𝑓 (𝑇𝑓𝑓 − 𝑇𝑙 ) 𝐿 − = 𝑉𝐻 𝐻 𝑎 𝜈𝑖 𝑘𝑓𝑓 (𝑇𝑓𝑓 − 𝑇𝑙 ) 𝑘𝑢𝑓 (𝑇𝐵𝑂𝑇 − 𝑇𝑓 ) 𝑑𝑧 − = 𝜌𝑠𝑖 𝐿 𝑎 𝑍 𝑑𝑡 Where a = thickness of the frozen fringe (m); H = thickness of frozen zone (m); kf = thermal conductivity of the frozen zone (W/(oC.m)); 38 Equation 2.51 Equation 2.52 kff = thermal conductivity of the partially frozen zone (W/(oC.m)); kuf = thermal conductivity of the unfrozen zone (W/(oC.m)); L = latent heat of fusion of water (J/Kg); TTOP, TBOT = temperatures at the top and bottom of the soil column (oC); VH = frost heave rate (m/s); Z = distance between bottom of soil column and position of ice penetration (m); Tl = temperature at the base of the active ice lens (oC); Tff = temperature at the base of the frozen fringe (oC); 𝑑𝑧 𝑑𝑡 = frost depth propagation rate (m/s); and 𝜌𝑠𝑖 = mass of ice per unit volume of soil (kg/m3). The water pressure at the Tf boundary is: 𝑃𝑤𝑓 = −𝑔 𝑃𝑤𝑓 𝑉𝑢𝑓 𝑍 (1 + ) 𝜈𝑖 𝐾𝑢𝑓 𝑑𝑧 𝑧 𝜈𝑤 (𝑉𝐻 + 𝜌𝑠𝑖 Δ𝜈 𝑑𝑡 ) = −𝑔 (1 + ) 𝜈𝑤 𝜈𝑖 𝐾𝑢𝑓 Equation 2.53 Equation 2.54 And finally the velocity in the active frozen zone can be obtained as: 𝑉𝐻 = 𝜈𝑖2 𝑔𝜈𝑤 𝐿(−𝑇𝑙 ) [ − 𝑃𝑂𝐵 + 𝑃𝐿𝑓 ] 𝑎𝐼𝑓𝑙 𝜈𝑤 𝑇𝑚 1 [𝑇 − 𝑇 ] + (𝐾 ) 𝑓𝑓 𝑙 𝐿 1 Equation 2.55 𝑇𝑙 1 𝑑𝑇 𝑇𝑓𝑓 𝐾𝑓𝑓 𝐼𝑓𝑙 = ∫ Where 𝑔 = acceleration of gravity (m/s2); 𝑃𝑤𝑓 = water pressure at the edge of the frozen fringe (kPa); Δ𝜈 = specific volume difference (νi − 𝜈𝑤 ); POB = overburden pressure (kPa); 39 Equation 2.56 KL = permeability of ice lenses (m/s); 𝐾𝑢𝑓 = permeability of unfrozen zone (m/s); and All other parameters are the same as before. Using different boundary conditions and solving Equation 2.51, Equation 2.52, Equation 2.54 and Equation 2.55 simultaneously, the water pressure at the bottom of the frozen fringe and the heave rate can be obtained. In addition, Gilpin proposed an approximate analytical solution. Nixon (Nixon, 1991) modified the approximate analytical solution of Gilpin. In this approach, a relationship between the frozen hydraulic conductivity and temperature is needed to predict the distinct location of each ice lens within the frozen zone. As shown in Figure 2.9, a linear temperature distribution (see Equation 2.57) and permeability distribution (see Equation 2.58) across the frozen fringe were assumed. 𝑥 𝑇 = (𝑇𝑙 − 𝑇𝑓𝑓 ) (1 − ) + 𝑇𝑓𝑓 𝑎 𝑘= Equation 2.57 𝑘𝑢𝑓 𝛼 𝑥 [−(𝑇𝑙 − 𝑇𝑓𝑓 ) (1 − 𝑎 ) − 𝑇𝑓𝑓 ] Equation 2.58 Where x = the depth from the face of the active ice lens (cm); All other parameters are as before. If the assumptions of no pore-water phase expansion and incompressible soil are made, the continuity of water flow indicates that 𝑑 =0 𝑑𝑃𝑤 𝑑𝑥 { } 𝑑𝑥 Equation 2.59 So the velocity of water flow should be constant in the frozen fringe. At any temperature the unfrozen water content can be characterized by 40 𝑤𝑢 𝐴(−𝑇)𝐵 𝑊𝑢 = = 𝑤𝑡𝑜𝑡 𝑤𝑡𝑜𝑡 Equation 2.60 Where T = the temperature (oC); 𝑤𝑢 = the gravimetric unfrozen water content; 𝑊𝑢 = the fraction of the unfrozen water content; 𝑤𝑡𝑜𝑡 = the total gravimetric moisture content; and A , B = constants. Figure 2.9 a zone of frozen soil, a freezing fringe, and an underlying zone of unfrozen soil (Nixon ,1991) The frost heave can then be calculated as 41 𝑎 𝐻𝑓 = 𝑛 ∫ (1 − 𝑊𝑢 )𝑑𝑥 Equation 2.61 0 Where n = porosity of soil; Hf = frost heave (cm); and All other parameters are as before. The unfrozen water content parameter is redefined as follow: Al = A/𝑤𝑡𝑜𝑡 , and the distribution of 𝑊𝑢 with depth in the frozen fringe x is 𝑊𝑢 = 𝐴(−𝑇)𝐵 = 𝐴[−(𝑇𝑙 − 𝑇𝑓𝑓 )(1 − 𝑥/𝑎) − 𝑇𝑓 ]𝐵 Equation 2.62 After integration, the frost heave can be calculated using Equation 2.63 1+𝐵 1 + 𝐴{(−𝑇𝑓 ) − (−𝑇𝑙 )1+𝐵 } 𝐻𝑓 = 𝑛 𝑎 (1 + 𝐵)(𝑇𝑓𝑓 − 𝑇𝑙 ) Equation 2.63 The frost heave rate can be obtained from Equation 2.64 𝑑𝐻𝑓 𝐿𝑑𝐻𝑓 𝑑𝐻𝑓 𝑑𝑎 𝑑𝑇𝑙 𝑑𝑇𝑙 = 𝐿( 𝑑𝑎 + ) = (𝑄𝑓𝑓 − 𝑄𝑢 ) 𝑑𝑡 𝑑𝑎 𝑑𝑡 𝑑𝑡 Equation 2.64 Where Qff = heat flux through frozen fringe; Qu = heat flux through unfrozen zones; and All other parameters are the same. By comparison with a numerical solution, the assumptions of linearity of the temperature profile can be checked. The results of the finite difference calculation and the comparison with the approximate analytical solution are displayed in Figure 2.10. Also the model was used for different kinds of soil and the comparison between the predicted and observed laboratory results was made as shown in Figure 2.11 for one of the cases. 42 Figure 2.10 Numerical verification of linear temperature profile assumption. T = + 2.7 oC at sample base, cooling rate = 0.84 oC /day, sample height = 10 cm, initial lens temperature = 0.14"C, A = 0.05, B = 0.5, W = 20%, and initial freezing point = - 0.04 oC, (Nixon 1991) Figure 2.11 Predicted and observed heave for Konrad test No. 4 (Nixon, 1991) 43 Fowler and Krantz (Fowler and Krantz ,1994) developed a generalized model for the secondary frost heave. In order to simplify the governing model’s equations, dimensional analysis techniques (i.e. normalization and scaling) were used and a new dimensionless parameter was introduced. The model could predict the thickness of ice lenses. It was also shown that the thickness of the frozen fringe initially starts to increase then it reaches a steady state. The results were in good agreement with experimental data based on a step freezing process. The model is also capable of being extended to incorporate the solute effects on the freezing temperature and the unsaturated soils effects in the secondary frost heave. Konrad and Morgenstern (Konrad and Morgenstern, 1980) proposed a semi empirical model to solve the mass and heat transfer. The model is based on two assumptions. The first is zero overburden pressure, which implies that the weight of the overlying soil can be ignored and the second, Clapeyron equation at the base of the ice lens is valid. 𝑃𝑤 = In term of total head; Neglecting the elevation head yields; 𝜌𝑤 𝐿 (𝑇 − 𝑇𝑚 ) = 𝑀𝑇𝑙∗ 𝑇𝑚 𝑃𝑤 𝐻𝑤 = ( ) + ℎ𝑒 𝛾𝑤 𝐻𝑤 = ( 𝑃𝑤 𝑀 ) = ( ) 𝑇𝑙 𝛾𝑤 𝛾𝑤 Equation 2.65 Equation 2.66 Equation 2.67 Where 𝑇𝑙∗ = 𝑇 − 𝑇𝑚 (oC); 𝑀= constant; ℎ𝑒 = elevation head (cm); 𝐻𝑤 = total head (cm); Tl = the temperature at the bottom of the frozen zone (the top of the frozen fringe) which depends on the soil type (oC); and 44 All other parameters are as before. Assuming Darcy law is valid and considering a two layered system consisting of unfrozen soil of thickness lu(t) having hydraulic conductivity ku and a frozen fringe thickness d(t) with the overall hydraulic conductivity ̅̅̅ 𝐾𝑓 (𝑡), the velocity of water movement can be attained using Equation 2.68. 𝑣(𝑡) = |𝐻𝑤 (𝑡)| 𝑑(𝑡) 𝑙𝑢 + ̅̅̅ 𝐾𝑢 𝐾𝑓 (𝑡) Equation 2.68 Where lu(t) = thickness of unfrozen soil (cm); ku = hydraulic conductivity in unfrozen zone (cm/s); ̅̅̅ 𝐾𝑓 (𝑡) = overall hydraulic conductivity of frozen fringe (cm/s); d(t) = frozen fringe thickness (cm); 𝑣(𝑡)= water flow velocity (cm/s); and All other parameters are the same as before. At last, the integration of the heave rate over the duration of freezing yields the total heave hs(t) as stated in Equation 2.69. 𝑡 ℎ𝑠 (𝑡) = ∫ 0 𝑡 𝑑ℎ𝑠 𝑑𝑡 = 1.09 ∫ 𝑣(𝑡)𝑑𝑡 𝑑𝑡 0 Equation 2.69 For one-dimensional heat flow the above Fourier equation can be written as 𝜕 𝑑𝑇 𝜕𝑇 (𝑘 ) + 𝑄 = 𝐶 𝜕𝑧 𝑑𝑧 𝜕𝑡 Equation 2.70 where C = volumetric heat capacity (J/(kg.oC)); k = thermal heat conductivity (W/(m.oC)); Q = internal heat generation term per unit area and per unit time (J/(s.kg)); and All other parameters are the same as before. 45 By solving the coupled heat and mass flow equations, heave can be calculated. The 𝑇𝑙 and ̅̅̅ 𝐾𝑓 (𝑡) are physical parameters of the soil that can be determined in the laboratory and must be known in order to use equations Equation 2.65 to Equation 2.70. In another development, Konrad and Morgenstern introduced the concept of segregation potential (SP). They conducted a simple linear analysis based on the following three assumptions 1. Clapeyron equation is valid at the ice lens base. 2. Water flows continuously with an overall hydraulic conductivity in frozen fringe. 3. The temperature (Tl) at the top of the frozen fringe measured in the laboratory for certain soil type is the same as that in the field. This temperature is called segregation temperature. The results of their analysis indicate that when the temperature of the warm-side of a soil sample is held constant and the other side freezes under various freezing temperatures, the water intake flow toward the ice lenses (heave rate) increases linearly as the temperature gradient increases. The slope of such linear line is called segregation potential (SP) which can be expressed as: 𝑉𝐻 = 𝑆𝑃 ∆𝑇 Equation 2.71 Where ∆T = the temperature gradient in the frozen fringe (oC/cm); SP = segregation potential (cm2/(day. oC)); and all other parameters are the same as before. They also conducted laboratory test in order to evaluate the theory. Their laboratory results were consistent with the theory (Konrad and Morgenstern, 1981). Nixon 1991 stated that the segregation potential theory published by Konrad and Morgenstern (Konrad and Morgenstern, 1981) addressed the velocity of the migrating water toward the freezing front and temperature gradient in the frozen fringe. He stated that “the 46 velocity of water arriving at an advancing frost front is related to the temperature gradient in the frozen soil just behind the frost front”. Since this theory is empirical, laboratory tests must be conducted to find the segregation potential in each soil type under different field conditions. Also, since laboratory and field conditions are not necessarily based on the same physical conditions, the predictions might not be fully reliable (Gilpin, 1982). Additionally, Konrad and Morgenstern investigated the effect of the overburden pressure on the frost heave rate. They concluded that 1. As the overburden pressure increases, the segregation temperature decreases. 2. Increasing the overburden pressure causes decrease in the unfrozen water content and consequently decrease in permeability. They reported that 400kPa overburden pressure causes 25 percent decrease in the permeability relative to zero overburden pressure. 3. Decrease in segregation temperature leads to lager frozen fringe. 4. Increasing overburden pressure causes decrease in the heave rate. The reason is that increases in the overburden pressure cause decrease in the overall permeability and decrease in the suction pressure to move water toward the frozen front. They also, investigated the concept of “shut-off pressure” at which no water will flow into or out of the soil. This concept is controversial, some researchers showed that such pressure can be found in different soil in laboratory conditions and water is drawn to the freezing front in pressures less than the shut-off pressure and expelled from the frost front as the pressure exceeded the shut-off pressure. Others believe that given a sufficient freezing time, the water expulsion could be followed by water intake again. 47 Konrad and Morgenstern also pointed out that because frozen fringe is relatively large in the field and frost penetration rate is small it is reasonable to assume that the Tf/dt is zero. Applying this condition in laboratory, they developed an equation for SP and for different overburden pressure (POB) as follows 𝑆𝑃 = 𝑎 exp(−𝑏𝑃𝑂𝐵 ) Equation 2.72 Where a and b = statistical constants that can be obtained by modeling the data obtained from laboratory tests; and All other parameters are the same as before. Gilpin (Gilpin, 1982) tried to relate the SP approach to his model (Gilpin, 1980). This model gave a physical foundation to the empirical model. He introduced a dimensionless segregation potential (DSP) as follows: 𝐷𝑆𝑃 = 𝐿 𝑉𝐻 𝜗𝑖 𝐾𝑓𝑓 (𝑇𝑙 − 𝑇𝑓 ) Equation 2.73 Where DSP = dimensionless segregation potential; and All other parameters are the same as before Then, he used Equation 2.73 and the laboratory data to obtain the constant values in the following two equation forms 𝐷𝑆𝑃 = 𝐶1 (𝑃𝑂𝐵 + 𝑃𝑠𝑒𝑝 − 𝑃𝐿𝑓 )−𝛼 Equation 2.74 −(𝑃𝑂𝐵 − 𝑃𝐿𝑓 ) 𝐷𝑆𝑃 = 𝐶2 exp [ ] 𝑃𝑠𝑒𝑝 Equation 2.75 Where all parameters are the same as before. The results indicated that the correlation coefficients are approximately the same (0.97) for both forms, so either one of them could be used for the calculation of DSP. 48 The applicability of segregation potential was investigated by some researchers (Nixon, 1982; Hayhoe and Balchin, 1990). Nixon installed two circular frost heave test plates at Foothills Pipe Lines test facility in Calgary, Canada and compared the measured data to the calculated ones. The SP values were obtained from laboratory data. The water intake at the bottom of ice lens was calculated using the following equation 𝑉𝑓𝑓 = 𝑆𝑃 𝑔𝑟𝑎𝑑 𝑇 Equation 2.76 Where Vff = the velocity of water flow in the frozen fringe (cm/day); SP = segregation potential (cm2/(day.oC)) and; grad T= the temperature gradient in the frozen fringe (oC/cm). The total heave was then estimated the following three equations ∆ℎ𝑡𝑜𝑡𝑎𝑙 = ∆ℎ𝑠 + ∆ℎ𝑖 Equation 2.77 ∆ℎ𝑠 = 1.09 ∗ 𝑉𝑓𝑓 ∗ ∆𝑡 Equation 2.78 ∆ℎ𝑖 = 0.09 ∗ 𝑛 ∗ 𝐻 Equation 2.79 Where ∆ℎ𝑡𝑜𝑡𝑎𝑙 = total frost heave (cm); ∆ℎ𝑠 = frost heave due to water intake (cm); ∆ℎ𝑠 = frost heave due to in-situ pore water freezing (cm); ∆𝑡 = time interval (day); N = soil porosity and; and H = frost depth (cm). The results were in relatively good agreement with the measured data in the field (Nixon, 1982). Other researchers used field data in Ottawa, Canada and the segregation potential approach to calculate the frost heave (Hayhoe and Balchin, 1990). Their data showed that SP 49 could change due to the change in soil properties with depth. However assuming a fixed value for SP, the estimated heave values were in relatively good agreements with the field data. Han and Goodings (Han and Goodings, 2006) investigated the differences between clay and silt freezing behavior due to lower hydraulic conductivity and higher water content of the saturated clay. They used a geotechnical centrifuge in order to observe the behavior of the soil. They found that in the unfrozen zone, consolidation occurs due to water migration toward the frozen zone. This consolidation reduced the total heave. The results of their tests indicated that, as for other soil type, the heave in clay decreases with increasing overburden pressure. They also found that due to the low hydraulic conductivity of the clay, the water content effect heave more than the ground water level (GWL). In other words, due to low hydraulic conductivity, it requires long time for water to migrate from the GWL toward the frozen zone. Therefore freezing in clay appears to be a close system and the immediate accessible supply of water has more effect than the GWL. Further, they developed a simple analytical model (see Equation 2.80) based on using consolidation concept for 100% saturated soils. 0.09𝑒𝑓 ∆𝐻 = ( ) ∗ 𝐻𝑓 1 + 1.09𝑒𝑓 Equation 2.80 Where ∆𝐻= heave (cm); 𝐻𝑓 = frost depth including heave (cm); and 𝑒𝑓 = final void ratio in the frozen fringe. For the case where saturation falls below 100% during the freezing the heave can be estimated as follows 1.09𝑒𝑓 − 𝑆𝑟𝑓 𝑒𝑖 𝑒𝑓 − 𝑒𝑖 ∆𝐻 = [( )−( )] ∗ 𝐻𝑓 𝑆𝑟𝑓 + 1.09𝑒𝑓 1 + 1.09𝑒𝑓 Where 𝑆𝑟𝑓 = degree of saturation in the frozen zone, 50 Equation 2.81 𝑒𝑖 = initial void ratio. Frost Pressure The magnitudes of heave pressure and heave caused by frost action vary from one scenario to another and they are function of water availability, soil type, overburden pressure, below freezing temperatures and the duration of cold season (Loch and Miller, 1972; Nixon, 1991; Konrad and Morgenstern, 1981). The magnitude of frost pressure highly depends on the particle sizes. In fine sand, the pressure is low, whereas it is intermediate for silt and high for clay. In fact, the pressure could vary between 420 psf in sand to 6300 psf in clay (Hoekstra, 1969). For certain design scenario where the frost potential cannot be eliminated, the heave pressure and heave should be accounted for in the design phase of the structure. Therefore, accurate estimation of the magnitudes of frost pressure and frost heave are essential part of the design process. The direction of frost heave and frost pressure is parallel to the heat flow direction (Penner and Irwin, 1969). For example, under the pavement and bridge foundations, the frozen front advances vertically downward whereas for retaining wall, frost heave progresses vertically and horizontally behind the wall (Andersland and Anderson, 1978). In order to eliminate or decrease frost heave potential, building foundations are typically placed below the frost line. Otherwise, frost will cause upward pressure against the foundation. When the upward pressure becomes higher than the downward pressure (due to the foundation weight and applied load), the foundation will move upward as shown in Figure 2.12a. Behind most retaining structures, the frost pressure is oriented horizontally against the wall as shown in Figure 2.12b. The combination of frost and active earth pressures may cause the wall to slide horizontally along its foundation. 51 The behavior of frozen soil was of interest in the past century, particularly after the ground freezing techniques were developed. Vialov (Vialove 1965) was the first one who investigated the viscoelastic behavior of frozen soil comprehensively. The viscoelastic behavior is especially important in the estimation of creep in artificial ground freezing (Lackner et al., 2008; Klein, 1981). Estimation of frost pressures and relaxation is of interest in cold-regions tunnel design. Klein (1981) developed a finite element time dependent model for investigating the behavior of temporary frozen earth support system in tunneling. Lai et al. (2000) used numerical inversion of Laplace transform for calculating forces and lining stress in tunnels. They assumed elastic-viscoelastic behavior for the frozen rock and used Poyting-Thomson model for illustrating the viscoelastic behavior in their model. Yuanming et al. (2005) assumed that the frozen soil is nonlinear elastic- plastic isotropic body and developed a finite element model for estimating the frost heaving pressure along a pile of a land bridge in china. Unfortunately, none of the mentioned researches conducted laboratory or field testing to evaluate their results. In the field, there are different complications in estimating frost pressure. The pressure could be different when the amount of free water is different or when the level of homogeneity is different. Sometimes the ice penetrates along the cracks and fissure instead of through the pores which results in lower frost pressures (Penner and Irwin, 1969). Different researchers reported measurements of frost pressure in different conditions and for different soil types (Penner 1969; Penner and Gold 1971; Kinosita 1967). Kinosita measured frost pressure both in field and laboratory conditions. His results showed decrease in frost pressure when the frost depth penetration stops or slows down in the cycles of freeze and thaw. He assumed that the decrease is because of viscoelastic behavior of frozen soil; i.e. stress relaxation. 52 (a) (b) Figure 2.12 Frost action under foundation causing uplift pressure and behind retaining structure causing horizontal pressure The results of field testing showed the following equation for decrease in frost force after a stop in frost penetration 𝐹 = 𝐹0 (𝑎𝑡 + 1)−𝑛 53 Equation 2.82 Where F0 = the force at the time of stoppage (Pa); t = time from frost penetration stoppage (hour); a and n = constants that can be obtained from experiments ; for temperature of -4 oC, a and n are 35 hour-1 and 1/6, respectively. By assuming a viscoelastic relationship between frost heave rate and force, he also suggested the following equation for estimation of frost heave force 𝐹(𝑡) = 𝐶𝑏 𝑎𝑡(2𝑛 − 1) + 1 {1 − } (𝑎𝑡 + 1)𝑛 𝑎(1 − 𝑛) Equation 2.83 Where b = constant frost heave rate; C = constant; and All other parameters are the same as before. His laboratory results also showed that frost heave and frost pressure have very similar trend. 2.3 Thaw Depth and Seasonal Load Restrictions This section consists of review of the literature regarding the thaw depth during spring season and the load restriction models. Thaw Depth With the assumptions of no heat transfer in the frozen zone and a linear temperature distribution in the thawed zone, to estimate the thaw depths over time, Stefan simplified Neumann’s equation as follow (Jiji, 2009) 𝑋=√ 2𝑘𝑢 𝑇𝑠 𝑡 𝜌𝑙 Where X = thaw depth (cm); t = time since the thawing starts (s); 54 Equation 2.84 ku = thermal conductivity of the unfrozen soil (W/(m.oC)); 𝑇𝑠 = applied constant surface temperature (oC); 𝑙 = latent heat of fusion (J/Kg); and 𝜌 = density (Kg/m3). Nixon and McRoberts (Nixon and McRoberts, 1973) modified the Stefan solution for multi layered systems. They considered a two layered system, the first layer with the depth of H and thermal conductivity of k1 and volumetric latent heat of L1 overlying the second layer of soil with thermal conductivity of k2 and volumetric latent heat of L2. They calculated the time to thaw the upper layer completely using Equation 2.85. Re-arranging Equation 2.85 yields Equation 2.86. The product 𝑡0 𝑇𝑠 in the last equation is defined as the cumulative thawing degree day (CTDD) that is required to completely thaw the first layer. Further, they estimated the thaw depth for the second layer using Equation 2.87 𝐻 2 𝐿1 2𝑘1 𝑇𝑠 𝑡0 = 𝑡0 𝑇𝑠 = 𝐶𝑇𝐷𝐷 = 𝑋 = [( Equation 2.85 𝐻 2 𝐿1 2𝑘1 𝑘2 2 𝑘2 𝑘2 𝐻) + 2 𝑇𝑠(𝑡 − 𝑡0 )]1/2 − ( − 1) 𝐻 𝑘1 𝐿2 𝑘1 Where the subscript 1 refers to the first layer; t0 = time to thaw the overlying layer (s); H = thickness of the first layer (m); L = volumetric latent heat of fusion (J/m3); k = thermal conductivity (W/(m.oC)); Ts = the mean surface temperature during the thawing period (C); 55 Equation 2.86 Equation 2.87 X= the thaw penetration depth (m); the subscripts 1 and 2 refer to first and second layer, respectively; and All other parameters are the same as before. Thus, the estimation of thaw depths using Equation 2.85 through Equation 2.87 requires knowledge of the thermal conductivity and latent heat of fusion of the soil. Since these inputs are not always available or expensive to collect, average values are typically used in the calculations of frost depths. Seasonal Load Restriction A typical pavement structure consists of two or three layer system depending on the pavement class. The thicknesses of these layers vary substantially from one road class to another. For all Interstate and primary roads, the thicknesses of the layers are designed and constructed to provide adequate protection of the roadbed soil against freezing. Such protection is not provided for the majority of the secondary roads including most county roads, city streets and farm-tomarket roads. During the winter season, available water in the pavement structure freezes creating ice lenses and causing increases in the stiffness of the various pavement layers. Over the winter months, the ice lenses grow in volume due to migration of water from the ground water table toward the freezing front. The growth of ice lenses causes the pavement to heave. During spring season, the frozen ice lenses melt from the top due to warmer temperature and from the bottom due to the internal heat of the earth. Melted water at the top of the ice lenses cannot drain by gravity because of the impermeable ice. Hence, the water from the melted ice saturate the pavement layers especially the upper portion of the roadbed soil. This causes substantial softening of the roadbed soil. Figure 2.13 shows the pavement deformations due to frost heave and thaw consolidation. At the beginning of the winter season, suction is build up at the frost 56 front due to frost action. This suction can cause volume reduction in the subgrade soil (See path “a”). On the other hand, as the frost front progresses, the suction causes upward water flow to the frost front which leads to formation of excessive ice lenses. As a result of ice lenses formation, the subgrade soil volume increases as shown in path “b”. When spring thaw occurs, the ice lenses start to melt from top down making the subgrade soil saturated and resulting in bearing capacity reduction (See path “c”). The extent of effective stress reduction depends on the degree of saturation, rate of thaw and drainage system efficiency. On the other hand, as thaw begin, excess pore pressure start to dissipate gradually resulting in thaw consolidation. Also, as the water drains out, depends on the subgrade soil permeability and effective hydraulic gradient the bearing capacity recovers slowly (See path “d”). It is also possible that freeze and thaw action results in negative of positive residual volume change at the end (Dore, 2004). ∆h h ∆𝝈’ 𝝈’ Figure 2.13 Frost heave and thaw consolidation process (Dore, 2004) 57 It is at this critical time, the damage delivered by the traffic load increases substantially causing the well- known spring break-up of the pavement structures. Various studies have shown that up to ninety percent of the pavement damage occurs during the yearly spring thaw period (Tighe et.al, 2007; Janoo, 2002). Table 2.3 depicts the changes in the MR and the resulting damage in a typical year (AASHTO, 1993). It can be seen from the table that the resilient modulus of the roadbed soils decreases during the thawing period and the relative damage increases substantially. Table 2.3 Yearly variation in the resilient modulus of roadbed soils and the associated damage (AASHTO, 1993) Month Roadbed Soil Relative Modulus, MR Damage (Kpa) uf January 137895 0.01 February 137895 0.01 March 17237 1.51 April 27579 0.51 May 48263 0.13 June 48263 0.13 July 48263 0.13 August 48263 0.13 September 48263 0.13 October 48263 0.13 November 27579 0.51 Decemeber 137895 0.01 To minimize this damage and to extend the life of the pavement structures, most highway authorities post seasonal load restriction (SLR) signs. A study in the State of Minnesota showed that 20% and 50% weight reduction are expected to increase the pavement service life by 62% and 95%, respectively (Huen et.al, 2006). The problem stems from the fact that the time for posting SLR varies substantially from one year to the next and it is a function of the environment 58 and the water flow regime in the area. Such practice is exercised by most European Countries, Canadian Provinces, and road owners in USA. Some agencies restrict the maximum load that can be carried by certain axles, others restrict load and speed. Since the severity of winter varies substantially from one year to the next, most existing practices unintentionally lead to some pavement damage. On the other hand, The SLR causes hardship to the trucking industry and increases the number of trucks on the road. Thus, accurate knowledge of the time when the SLR signs should be posted and removed is crucial to the road owners and the road users. Such timely posting and removing the SLR signs cannot be had unless accurate prediction of frost and thaw depths as a function of time can be accomplished. In addition, accurate prediction of freeze-thaw cycles is critical to an effective load restriction approach (Ovik et al., 2000). Various practices are used by road owners to establish the value of the maximum loads and the dates for posting and removing the SLR signs. These practices include: 1. Engineering judgment and visual observations. 2. Fixed dates for posting and removing the SLR signs. 3. Upper bound or critical values of pavement deflections measured using either a falling weight deflectometer (FWD) or a portable falling weight deflectometer (PFWD). Such practice is expensive and time consuming and requires repeated FWD testing (Kestler et.al, 2007, Tighe et.al, 2007). 4. Quantitative approaches based on location and observed severity of the winter season (Ovik et.al, 2000; Miller et.al, 2013, Kestler et.al, 2007, Tighe et.al, 2007). Kestler et.al (Kestler et.al, 2007) recommended using practice 3, 4 or combination of both along with complementary guidance stated in Table 2.4. Washington Department of Transportation (WSDOT) proposed a SLR guideline, which has been adapted and modified in 59 several State Highway Agencies (SHAs). In this policy, the CTDD is calculated using Equation 2.89 (WisDOT, 2003) 𝑛 𝐶𝑇𝐷𝐷𝑛 = ∑ 𝑇ℎ𝑎𝑤𝑖𝑛𝑔 𝑑𝑒𝑔𝑟𝑒𝑒 𝑑𝑎𝑦 Equation 2.88 𝑖=1 𝑇𝑚𝑎𝑥 +𝑇𝑚𝑖𝑛 𝑇ℎ𝑎𝑤𝑖𝑛𝑔 𝑑𝑒𝑔𝑟𝑒𝑒 𝑑𝑎𝑦 = ( 2 − 𝑇𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 ) Equation 2.89 Where Tmax = maximum daily air temperature (°C); Tmin = minimum daily air temperature (°C); and Treference = -1.7 °C. Since various temperature measurements indicated that during the thawing period, the asphalt pavement surface temperature is 0°C when the air temperature is about -1.7°C, the reference temperature was set at -1.7°C (Mahoney et.al., 1986). Table 2.4 provides a list of Washington’s seasonal load restriction policy for thin and thick pavements. The posting of the SLR signs are based on the two levels listed in Table 2.4; should and must is placed when the 5-day weather forecast shows that the CTDD will reach the “should be posted or the must be posted” levels. The CTDD value for each level and pavement type is listed in Table 2.4 (Mahoney et.al., 1986). Table 2.4 CTDD threshold for posting SLR (Mahoney et.al., 1986) Pavement Structure CTDD “Should” Level “Must” Level THIN Asphalt 2” or less Base course 6” or less 6°C- degree day 22°C-degree days THICK Asphalt more than 2” Base course more than 6" 15°C- degree days 28°C- degree days 60 They also developed a regression equation between freezing index (FI), CTDD and thaw duration (D) based on the results of a heat flow simulation model on fine-grained subgrades. They recommended to remove the SLR as soon as one of these conditions were reached (Mahoney et.al., 1986) 𝐷 = 25 + 0.006 ∗ 𝐹𝐼 Equation 2.90 𝐶𝑇𝐷𝐷 = 0.17 ∗ (𝐹𝐼) Equation 2.91 Where D = thaw duration (days); FI = freezing Index oC-days; and All other parameters are the same as before. The proposed thaw duration was consisting of the thaw duration and the base recovery time after thaw ends. Yesiller et al (Yesiller et al, 1996) showed that the WSDOT equations predicts the thaw duration for fine-grained better than granular materials. A study in Minnesota indicated that WSDOT equation does not accurately predict the thaw duration. Since WSDOT used numerical model for developing the equation, the model may not reflect the field conditions precisely and this might be a possible reason for deviation of the prediction from the actual duration (Ovik et.al, 2000). The Minnesota Department of Transportation (MnDOT) developed SLR policy based on the following criteria (Ovik et.al, 2000) 1. Thaw had penetrated to the depth of 15 cm. 2. The weather forecast shows continuous thaw. 3. Different district within the same zone should reach to the same thaw depth (Minnesota was divided into 5 zones). 61 It was also recommended that deflection measurement should be collected and SLR should be removed three weeks after the maximum deflection. The deflection measurements showed that maximum deflection occurred when the thaw depth was about 102 centimeters. Table 2.5 Complementary guild for SLR Which pavements need SLR? 1- Pavements with frost susceptible soils in base and subgrade. 2- Pavements with fine materials in subgrade such as ML,MH,CL and CH 3- Pavements with poor drainage system or high groundwater table. 4- Pavements in which distresses like fatigue cracking or rutting has been detected 5- Pavements in which surface deflection is 45% to 50% higher in spring than summer When SLR should be placed? 1- When pavement surface deflection increases up to 45-50% higher than summer. 2- When the thaw depths reaches the subgrade (This could be determined which a thaw index method or temperature measurements with frost tubes or electric resistance gauge) What load restrictions should be placed? 1- Allow load levels that reduce the pavement deflections to those occur in summer. 2- Allow load levels that lead to desired proliferation in service life. 3- Allow load levels based on guidelines developed by different researchers When SLR should be removed? 1-When pavements recover and the measured surface deflections reduce to the summer values. Henceforth, MNDOT revised the policy. They modified the thawing index calculation approach and adapted a variable reference temperature instead of a fixed one and also implemented SLR when thaw reaches 30.5 centimeters. The revised SLR policy is based on the cumulative thawing degree-days (CTDD) and the cumulative freezing degree days (CFDD) (MnDOT, 2009). The MnDOT procedure for the calculation of the CFDD and CTDD is summarized in a flowchart format shown in Figure 2.14. The decision to post and remove the SLR signs is based on the results of the calculation of the CFDD and CTDD and the corresponding reference temperatures are listed in Table 2.6. The reference temperature accounts for the effect of the duration and intensity of the sun 62 radiation on the pavement thawing. As stated before, the MnDOT’s SLR policy divides the state into five zones. The SLR is posted when the 3-day weather forecast shows that the CTDD for a given zone is more than 15 oC-degree day and the continued warmth is predicted for longer time period. CTDD0=0 and n starts at 1 For day n Taverage-T reference > 0 oC YES NO Thawing degree day = T average-T reference o CTDDn-1> (0.5*(0 C-T average)) NO Freezing degree day = 0 oC- day YES Thawing degree day = 0 oC- day Thawing degree day = 0 oC- day Freezing degree day = 0 oC- day Freezing degree day = 0 oC-T average n=n+1 CTDDn=CTDDn-1+(Thawing degree day- 0.5*Freezing degree day) Figure 2.14 MnDOT CDTT calculation flowchart. In this flowchart CTDDn is the cumulative thawing degree day calculated over ‘n’ days , CTDDn-1 is the cumulative thawing degree day calculated over ‘n-1’ days, Taverage is the daily average air temperature ((Tmax-Tmin)/2) and the Treference is the reference temperatures listed in Table 2.6 63 Table 2.6 Reference temperature Date Reference Temperature (oC) January 1 – January 31 0 February 1 – February 7 -1.5 February 8 – February 14 -2 February 15 – February 21 -2.5 February 22 – February 28 -3 March 1 – March 7 -3.5 March 8 – March 14 -4 March 15 – March 21 -4.5 March 22 – March 28 -5 March 29 – April 4 -5.5 April 5 – April 11 -6 April 12 – April 18 -6.5 April 19 – April 25 -7 April 26 – May 2 -7.5 May 3 – May 9 -8 May 10 – May 16 -8.5 May 17 – May 23 -9 May 24 – May 30 -9.5 June 1 – December 3 0 They also revised the thaw duration relationship as follow: 𝐷 = 0.15 + 0.010 ∗ 𝐹𝐼 − 19.1𝑃 − 12090 ∗ 𝑃 𝐹𝐼 Equation 2.92 Where P=frost depth (m); and All other parameters are the same as before. It is noteworthy that the standard error of estimate for this equation is 8 days and although the R2 is about 0.5, the equation predicts the thaw duration more accurately than WSDOT equations. However, because of the low correlation MnDOT, the length of SLR policy was considered to be 8 weeks. They also verified the 8 weeks length by using the rate of strength 64 recovery in back calculation moduli. It is of interest that their results showed that at two weeks past the end of thaw the pavement strength recovers between 50 to 100 percent depends on the soil type. The increase in fines in the base material leads to longer recovery period. Leong et al. (Leong et al, 2005) used Long Term Pavement Performance (LTTP) database in 6 stations in Ontario to develop a thaw index base method for predicting SLR. It should be noted that the frequency of the data in LTTP varies based on the data type, for instance the FWD and surface temperature data were collected once or twice every two years. They plotted the average daily temperature against the asphalt temperature and adapted the horizontal intercept of the graph as their reference temperature for the thaw index calculation. They analyzed the data for the first 100 oC-days of the year and obtained thawing index of 13 oC-days as the threshold. This threshold is very similar to MNDOT threshold thawing index of 15 oCdays. Berg et al (Berg et al, 2006) developed an alternative method for SLR, which considers the effect of pavement temperature. They assumed that the air temperature can be fitted into the following equation (See Figure 2.15) 2𝜋 𝑇𝑡 = 𝑀𝐴𝑇 + 𝐴𝑚𝑝 ∗ sin(( ) ∗ (𝑡 − 𝐿𝑎𝑔)) 𝑃 Equation 2.93 Where Tt = sinusoidal temperature on Julian day t; P = period of sinusoidal variation (365 days); MAT = mean annual temperature in 20 years; Amp= amplitude of temperature sinusoid and; Lag= the time that it takes for temperature sinusoid to reach the MAT (110 days Lag was used) 65 Surface Temp (oC) Air Temp (oC) Temperature (oC) Mean Annual Surface Temp (oC) Mean Annual Air Temp (oC) Freeze (oC) Julian Date (Days) Figure 2.15 Sinusoidal temperature variation (After Berg et al) They also estimated the surface temperature using the n-factor of 0.5 and 1.7, respectively. Then, they calculated the difference between air and surface temperature for each day. This difference was added to measured air temperature to estimate the surface temperature. The estimated surface temperature was used to calculate the surface thawing index (CTI) as follows 𝑛 𝐶𝑇𝐼 = ∑ 0 − (𝑝𝑎𝑣𝑒𝑚𝑒𝑛𝑡 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 𝑡𝑒𝑚𝑝𝑎𝑟𝑎𝑡𝑢𝑟𝑒) Equation 2.94 𝑖=1 They recommended to place the SLR with CTI reaches 17 oC-days but they did not provide a protocol for removing the SLR. It should be noted that this method does not work well in the mountain area where the closes weather station is not relatively close (Kestler et al, 2011). Bradley et al (Bradley et al. 2012) developed a SLR policy in Manitoba by using the FWD data, moisture content, and temperature profile in pavement layers. They calculated thawing index as follows 66 𝑇𝑚𝑎𝑥 + 𝑇𝑚𝑖𝑛 𝑇𝑚𝑎𝑥 + 𝑇𝑚𝑖𝑛 ≥ 0−→ 𝐶𝑇𝐼 = ∑(𝑇𝑟𝑒𝑓 + ) 2 2 Equation 2.95 𝑇𝑚𝑎𝑥 + 𝑇𝑚𝑖𝑛 𝑇𝑚𝑎𝑥 + 𝑇𝑚𝑖𝑛 < 0 −→ 𝐶𝑇𝐼 = ∑(𝑇𝑟𝑒𝑓 + ) 2 4 Equation 2.96 𝑖𝑓 𝑖𝑓 Where Tref = 1.7 oC at March 1 and increases by 0.055 oC per day until May 31 (CTI cannot be negative, i.e. it resets when CTI values drops below 0), and All other parameters are the same as before. They recommended a CTI threshold of 16 oC-days for posting the SLR. Based on their FWD data they found that the pavement recovers noticeably when thaw depth reaches to 119 cm. They suggested removing the SLR in 8 weeks or when CTI reaches 350 oC-days, whichever occurs first. Miller al (Miller et al, 2013) used Bradley et al, WSDOT and MnDOT methods and compared their results with measured data in New Hampshire. They used FWD data and temperature logger data in their analysis. Since Manitoba is closer in latitude to Minnesota, they calculated thawing index based on MnDOT method. Their results indicated that in posting SLR the MnDOT and Bradley et al methods are slightly conservative, whereas WSDOT method dates for applying SLR are consistently late. On the other hand, in most cases Equation 2.91 and Bradley et al threshold were able to predict the time of SLR removal accurately. However, both of these methods can be non-conservative in some cases. Generally speaking, the 8 weeks duration seems to yield fairly conservative results for SLR removal. Manitoba department of transportation developed the following Equation for thawing index calculation: 𝑇𝐴𝑖𝑟 2 𝑖𝑓 𝑇𝑀𝑒𝑎𝑛 < 0 𝑇𝑀𝑂𝐷 = 𝑖𝑓 𝑇𝑀𝑒𝑎𝑛 ≥ 0 𝑇𝑀𝑂𝐷 = 𝑇𝑀𝑒𝑎𝑛 67 Equation 2.97 Equation 2.98 𝑇𝐼𝑖 = 𝑇𝐼𝑖−1 + 𝑇𝑀𝑂𝐷 + 𝑇𝑅𝑒𝑓 Equation 2.99 𝑇𝐼𝑅𝑒𝑓 = 1.7 + 0.06𝑖 Equation 2.100 Where TMOD = Modified Temperature (oC); TAir = Air temperature (oC); TRef = Reference Air temperature (oC) – 1.7 (oC); TIi = Thaw index for day i (oC-day); TIi-1 = Thaw index for day i-1 (oC-day); and i = number of days since March 1. They recommended posting SLR when thawing index reaches 15 (oC-day); Kestler et al (Kestler et al, 2007) investigated the applicability and accuracy of three methods for posting and removing SLR and reached the following conclusions 1. Using subsurface temperature and moisture sensors: Their results showed that the temperature-moisture system can be reliable method to determine when to post and remove SLR signs. SLR can be commenced when the surface temperature nears 0 oC and can be ended when the moisture dissipates. It should be noted that although this method is fairly accurate, it is site specific. 2. Using the lightweight or potable FWD (PFWD): Since the FWD initial cost is beyond the some agencies budgets, the use of PFWD for placement and removal of SLR was investigated. Their results showed that the trend in composite moduli obtained from PFWD and the moduli back calculated using FWD data is the same. On the other hand, as the asphalt thickness decreases the correlation between the PFWD and PFW modulus increases. They concluded that PFWD can be used as a complementary approach for posting and removing SLR in low volume road with gravel or thin asphalt surface. 68 3. Using thaw index methods: the thawing index approach proposed by MnDOT showed promising results but required different reference temperatures that might vary by time and location. On the other hand, Since the Berg et al approach does not require a reference temperature; the method can be used in any site without calibration. However, their results indicated that although the method works well for posting the SLR it is unable to accurately predict the SLR removal. Furthermore, Kestler el at ( Kestler el at, 2011) developed a toolkit for implementing the SLR using the following three methods 1. Everseries suite: Everseries suite developed in the State of Washington was used to post SLR. Firstly, pavement moduli were back-calculated using Evercalc. Secondly, the stress and strain and damage factor were estimated using Everpave. Damage factor is the ratio of the load cycle required to reach failure in normal condition to the load cycle required to reach failure under thaw condition. They results showed the highest damage factor when thaw reaches 46 centimeters. Also, the damaged factor evened out after 5 weeks of implementing the SLR. 2. Enhanced Integrated Climatic Model (EICM): EICM was used as a tool to develop the thaw predictor. EICM was modified to predict the thermal and mass transfer between pavement layers. They indicated that since the output of the model is frost and thaw depth, therefore the model can be used as a tool to help the DOTs in the SLR decision making. They recommended using another climatic database along with National Climatic Data Center (NCDC). NCDC is consisting of 851 stations across the United States but most of these stations are at the airports or cities therefore there are not good representatives of remote site climate conditions. The recommended database is consisting of 2000 remote automated weather stations (RAWS) which are frequently at isolated areas. Also, since the low volume 69 road could be covered by few centimeters of snow, they also recommended considering the effect of snow cover in the model. Their results indicated that the EICM model can predict the beginning of the thaw season fairly well but still needs additional modifications to predict the end of the thaw as well. 3. The Lightweight Deflectometer (LWD): The LWD estimates the in-situ stiffness modulus of the soil layers. They showed that the LWD followed the seasonal trend in stiffness changes and has a good correlation with FWD-derived moduli. For posting the SLR, they recommended to take different tests in order to obtain modulus in the normal condition and then take multiple readings during spring thaw season. It was recommended to implement the SLR when composite modulus is less than 20% of the modulus in normal condition and end the SLR when the composite gain the 80% of modulus value. Chapin et al (Chapin et al., 2012) utilized finite element program TEMP/W (GEOSLOPE 2007) to simulate freezing and thawing front in the pavement. They applied the program to two sites in northern Ontario with considerably different pavement structures. First, a steady state analysis was conducted to establish the initial conditions within the model and second, a transient analysis was conducted. For the upper boundary conditions, they investigated two different scenarios. First scenario was converting the average daily air temperature to pavement surface temperature by using n-factor. They indicated that n-factor can be between 1.4 to 2.3 and between 0.29 to 1 for thaw and frost depth prediction, respectively. The second scenario involved using increasing weekly reference temperatures similar to MnDOT method. Their results showed that the first scenario yielded the most accurate predictions. But the n-factor can be different in different locations and even in different years. In Ontario n factor of 2.3 and 0.6 yielded the most accurate thaw and frost depth predictions, respectively. Their results indicated 70 that TEMP/W is able to accurately simulate the thaw depth and rate of thawing at the beginning of thaw, which is the critical time for posting SLR. But as the thaw progresses the accuracy of the simulation can be decreased. Therefore, they concluded that TEMP/W simulation results can be used along with the thawing index approach for posting the SLR but not for removing it. Moreover, they recommended that thawing index threshold for applying SLR should be refined actively as more data become available and also, considering the climate difference as well as pavement structures, different CTI threshold in different regions should be used. 71 CHAPTER 3 DATA MINING 3.1 Database The evaluation of the accuracy of existing frost depths and heave models or the development of new accurate and representative ones require field data that represent the environment and the various pavement structures (Tighe et.al, 2007). Fortunately, such data were available and obtained from the Michigan Department of Transportation (MDOT) and the Minnesota department of Transportation (MnDOT). The following databases and their sources were used in this study: 1. Road Weather Information System (RWIS) for frost depth and weather data measured in the state of Michigan. It should be noted that the RWIS subsurface sensors do not measure the frost depth directly, they measure the subsurface temperature. In the analyses, it was assumed that the ground water freezes at 0oC. 2. National Oceanic and Atmospheric Administration (NOAA) weather data in the states of Michigan and Minnesota. 3. Minnesota Department of Transportation (MnDOT) database for frost depth measured in the state of Minnesota. 4. Michigan Department of Transportation (MDOT) database for frost heave measured in the state of Michigan. 5. Michigan State University Enviro-weather (MSU-EW) for weather data in the state of Michigan. 72 3.2 Frost Depth Data 3.2.1 The State of Michigan RWIS uses different technologies that collect, transmit and publish weather and road condition information. The weather data is collected by the environmental sensor station (ESS). In these stations the sensors collect and transmit weather and pavement data (US DOT, 2002). In general RWIS may encompass: 1. Meteorological sensors for measuring atmospheric pressure, temperature, relative humidity, visibility, wind speed and direction, and precipitation (amount and type). 2. Pavement sensors for measuring pavement temperature and condition (wet, dry, snow), subsurface pavement temperature, the amount and type of deicing chemical used on the pavement surface. 3. Pavement temperature and weather condition forecast based on the site (Boselly et al., 1993). In this study, the RWIS database that was provided by MDOT was used for subsurface pavement temperature data (RWIS, 2012). RWIS consist of 25 stations located throughout the State of Michigan. However, in this study, only 18 stations were used (MDOT 2008, MDOT 2009a, and b) due to partially missing data in seven stations. Figure 3.1shows the stations’ location in the state of Michigan. Table 3.1 shows the RWIS stations ID, latitude, longitude and soil type. The detailed soil log for each station was provided by MDOT. In all stations, one year data (2010-2011) were available and used except for Au Train, Harvey and Brevort Stations in the UP where two years of data (2009-2010) were available and used. 73 Figure 3.1 RWIS station locations, Michigan As mentioned before the RWIS collect meteorological as well as surface and subsurface pavement temperature data. The data were collected approximately every 10 minutes. The RWIS data were provided by MDOT and processed in order to be used in this study. Figure 3.2 depicts a soil profile showing the locations of temperature sensors at a typical RWIS station. The data were then used to develop a GIS contour map of maximum frost depth in a typical year in Michigan as shown in Figure 3.3. 74 Table 3.1 RWIS stations in the State of Michigan Region Upper Peninsula Lower Peninsula Station Latitude Longitude Au Train 46.43 -86.84 Brevort 46.01 -85.01 Cooks Engadine 45.91 46.10 -86.48 -85.62 Golden Lake 46.16 -88.88 Harvey 46.49 -87.23 Michigamme 46.54 -88.13 Seney St. Ignace 46.35 45.90 -86.04 -84.74 Twin Lakes 46.88 -88.86 Benzonia 44.59 -86.10 Cadillac 44.25 -85.37 Grayling 43.89 -85.53 Houghton Lake 44.77 -85.40 Ludington 45.36 -85.18 Reed City 44.61 -84.71 Waters 44.33 -84.81 Williamsburg 45.76 -84.73 75 Subgrade Soil (up to 3.05 m) Upper 0.9 meter of Sand with Gravel and Silt, below Which Loose Moist Fine Sand Loose To Moderately Compact Moist Fine Sand Silty Clay with Sand Plastic Moist Sandy Clay Moderately Compact Moist Fine Sand with Gravel Loose to Moderately Compact Moist Fine to Medium Sand Upper 1.5 meter Clayey Sand, 1.5 meter and below Wet Find to Medium Sand with Silt Loose Moist to Wet Sand Silty Clay with Sand Loose Moist Fine to Medium Silty Clayey Sand Very Loose Fine to Medium Sand Medium Compact Fine Sand Fine to Coarse Sand with Gravel Loose to Medium Compact Fine Sand Very Loose Silty Sand with Trace Clay Loose to Medium Fine Sand with Trace Gravel Medium Compact to Loose Fine Sand with Trace Gravel Medium Compact Fine Sand with 0.5 meter of Silty Clay Unfortunately, the meteorological data were not available for the whole 2010-2011 winter. Therefore the NOAA and/or the MSU-EW weather data were used in the analyses (MSU-EW, 2012; NOAA, 2012). The selected NOAA and/or MSU-EW stations were within 16 kilometers of an RWIS station, otherwise the RWIS data were not used in the study. Table 3.2 shows the data availability in each database. Surface 0 cm 7.6 cm Base 13.2 cm 22.9 cm 38.1 cm 53.3 cm 71.1 cm Sub-base 86.4 cm 101.6 cm 116.8 cm 132.1 cm Roadbed Soil 147.3 cm 157.5 cm 172.7 cm 182.9 cm Figure 3.2 A soil profile at a typical RWIS station showing the depths of the temperature sensors 76 146.8 132.1 132.1 132.1 146.8 146.8 147.3 172.7 101.6 157.5 101.6 86.4 132.1 101.6 91.4 86.4 91.4 91.4 76.2 89.0- 93.0 cm 96.5- 101.6 cm 104.0- 106.7 cm 109.2- 111.8 cm 114.3- 119.4 cm 121.9- 124.5 cm 127.0- 129.5 cm 132.1- 137.2 cm 139.7- 182.9 cm Figure 3.3 Maximum frost depth contours in a typical year in the State of Michigan 77 Table 3.2 Data availability in each database Location Lower Peninsula Upper Peninsula 1. 2. 3. 4. RWIS Station Name Database Station Name Pr1 Temp2 WS3 SR4 Database Station Name Pr1 Temp2 WS3 Station Name Temp2 Pr1. Benzonia Benzonia ✓ ✓ ✓ ✓ N/A N/A N/A N/A N/A N/A N/A Cadillac McBain ✓ ✓ ✓ ✓ N/A N/A N/A N/A N/A N/A N/A Grayling N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A Houghton Lake Ludington MSUEW NOAA Ludington ✓ N/A ✓ Roben- Hood Airport Williamsburg Traverse City ✓ ✓ ✓ ✓ N/A Waters Gaylord ✓ ✓ ✓ ✓ Gaylord-9 SW Michigamme N/A N/A N/A N/A N/A N/A Harvey N/A N/A N/A N/A N/A Sawyer Int. Airport Au Train Chatham Brevort N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A Seney MSUEW McMillan ✓ ✓ ✓ N/A ✓ ✓ N/A ✓ N/A ✓ ✓ Reed City Engadine N/A ✓ Graying Army Airfied Roscommon Co. Airport ✓ ✓ N/A ✓ ✓ ✓ N/A ✓ N/A ✓ ✓ USC00203391 ✓ ✓ N/A ✓ ✓ N/A ✓ ✓ ✓ Big Rapids Waterwork ✓ ✓ N/A N/A N/A ✓ N/A N/A N/A ✓ N/A ✓ ✓ N/A ✓ N/A N/A N/A N/A N/A N/A N/A N/A N/A N/A Marquette ✓ N/A N/A N/A N/A N/A N/A N/A N/A N/A USC00205591 ✓ ✓ Engadine MDOT ✓ ✓ Naubinway ✓ N/A ✓ N/A ✓ N/A NOAA N/A Schoolcraft Co. Airport N/A ✓ N/A ✓ N/A Cooks N/A N/A N/A N/A N/A Twin Lakes N/A N/A N/A N/A N/A N/A ✓ N/A ✓ Stambaugh 2 SSE ✓ ✓ ✓ N/A Golden Lake N/A N/A N/A N/A N/A Kings Land O' Lakes Airport St. Ignace N/A N/A N/A N/A N/A St. Ignace Precipitation Temperature Wind Speed Solar Radiation 78 N/A N/A N/A N/A N/A Garden Corners N/A N/A ✓ N/A ✓ N/A N/A ✓ N/A ✓ N/A 3.2.2 The State of Minnesota In order to evaluate and verify the frost and thaw depth model which was developed using the measured frost depth data in the state of Michigan, the Minnesota frost depth data were requested, received and used in this study. The MnDOT data consisted of 9 years of data (2003 to 2012) collected at 8 stations located throughout the State of Minnesota as shown in Figure 3.4. Similar to the Michigan case, the meteorological database from the nearest NOAA station was obtained and used in the study. Figure 3.4 MNDOT stations location, Minnesota 3.3 Soil Properties Data Disturbed soil samples from different RWIS stations in the State of Michigan were provided by MDOT. The thermal properties of the soil were then measured in the laboratory at 79 Michigan State University. Since only disturbed soil samples were received, the insitu dry densities and water contents were unknown. Nevertheless, seven representative soil types were selected and their thermal properties were measured in the laboratory using KD2 pro thermal properties analyzer. The KD2 pro is a small and portable device with the capability of measuring different thermal properties of almost any material. The device has three sensors 6 cm single needle (KS-1), 10 cm single needle (TR-1), and 3 cm dual-needle (SH-1). Each sensor could be used depending on the thermal properties being measured and the material type. The KS-1 sensor was used to measure the thermal conductivity and thermal resistivity of the soil. The sensor is most accurate in liquid samples and insulating materials. In liquid samples free convection could be a source of error. Since the sensor applies small amount of heat to needles, free convection could be prevented making KS-1sensor a good choice for liquid samples. On the other hand, in granular samples like soil or powders, contact resistance could be a source of error. Size of the sensor and short heating time could maximize this error making the KS-1sensor a poor choice for these types of materials (Decagon Devices Inc., 2008). The TR-1 sensor was also used to measure the thermal conductivity and thermal resistivity. TR-1 is a hollow needle and comes with drill bits. The size and relatively longer heating time could minimize the contact resistance error, which make this sensor the primary choice for soil and granular materials. The sensor fully complies with ASTM D5334-08 specifications and standard test procedure for measuring the thermal conductivity of soils (Decagon Devices, Inc., 2008). The SH-1 sensor measures volumetric heat capacity, thermal diffusivity as well as thermal conductivity and thermal resistivity. This sensor could be used in most solid and 80 granular material but not in liquid samples. Table 3.3 summarizes the applicability of each sensor for different materials (Decagon Devices, Inc., 2008). Table 3.3 Sensor use guide (Decagon Devices Inc., 2008) Sample Material Low viscosity liquids (Water) High viscosity liquids (glycerol, oil) Insulation and insulating materials Moist soil Dry soil, powders, and granular materials Concrete and rock Other solids KS-1 Best Best Best N/A N/A N/A N/A TR-1 N/A OK N/A Best Best Best Best SH-1 N/A N/A N/A OK OK OK OK As stated before, seven types of soil were tested to obtain the required thermal properties. TR-1 and SH-1 sensors were used for measurements (Figure 3.5). In order to minimize the error five measurements were done on each sample as shown in Figure 3.6. The thermal properties of the samples were considered to be the average of these readings. Table 3.4 depicts the measured thermal properties of seven different types of soil in saturated condition. It is noteworthy that all soil samples were disturbed and were not compacted in the laboratory. Figure 3.5 Thermal conductivity measurement using KD2 pro 81 As mentioned before, detailed soil log for each RWIS station were available in the State of Michigan. On the other hand, the soils at the State of Minnesota stations were categorized as clayey and sandy soils. Therefore, in the model the average values of Table 3.4 for clayey and sandy soils were used for both states. 3.4 Frost Heave Data Data from 5 sites in Oakland County, Michigan were provided by MDOT (Novak, 1968). The data was collected in the winter of 1962-63 in a six-mile section (EBI 63172, CR5H) along I-75; Figure 3.7 depicts the location of the sites. Table 3.5 shows soil type, the maximum measured frost heave and frost depth under pavement and shoulder and at different stations. The detailed measurements with time can be found in Appendix B. Side of the mold Probe hole Figure 3.6 Five locations where the soil thermal properties were measured in each soil sample using KD2 pro 82 Table 3.4 Measured thermal properties of different types of soil using KD2 Pro Station Name Moisture Condition Material Williamsburg Silty Fine Sand with Trace of Gravel Fine Sand Fine Sand with Trace of Gravel Fine Sand Soft Clayey Sandy , Some Silt & Some Gravel Silty Clay Rudyard Silty Clay Houghton Lake Wolverine Saturated Thermal Conductivity (W/m.K) Heat Capacity (MJ/m3.K) 2.57 2.67 2.55 2.84 2.49 2.69 2.42 2.69 1.74 2.98 1.51 3.1 1.12 3.2 It should be noted that for station/528+88 the difference between the measured frost heave under the pavement and under the shoulder is approximately 1 centimeter which is much higher than the other sites. At this site, an undercut of approximately 30.5 centimeter was made while constructing the pavement for frost protection. 3.5 Pavement Surface Temperature and Thaw Depth Data In order to develop an effective spring load restriction policy, accurate models for calculating Cumulative Thawing Degree Day (CTDD) and estimating thaw depth are required. For developing a CTDD model, air temperature, pavement surface temperature data from the 18 RWIS stations in 2011, 2012, 2014, and 2015 were used (See Table 3.1). For developing a thaw depth model, data from 14 RWIS station in 2011, 2014, and 2015 were used. Unfortunately, the subsurface temperature data were not available for the whole 2012 spring. Therefore, that year was not included in the study. In addition, among the 18 RWIS stations in Table 3.1, the data from Cooks, St. Ignace, Twin Lake, and Williamsburg were excluded from the study due to partially missing data in 2011, 2014 and 2015. 83 Figure 3.7 MDOT frost heave station locations, Oakland County, Michigan. Table 3.5 Measured total heave and frost depths in different soil types, I75, Oakland County, Michigan Station Name Sta/724+00 Sta/719+00 Sta/652+00 Soil Type Fine Sand and Silt with Pebbles Fine Sand with Silt Pockets with Pebbles Insitu Sub Soil Clayey, Silty, Gravely, Sand Frost Depth (cm) Max Heave in Shoulder (cm) Max Heave in Pavement (cm) Duration (days) 60.96 2.54 1.91 65 71.12 2.16 1.91 40 86.36 2.29 2.16 60 Sta/528+88 Sandy Clayey Silt 76.20 1.91 1.02 70 Sta/474+00 Clayey Silt 63.50 2.54 2.29 55 84 CHAPTER 4 DATA ANALYSIS & DISCUSSION Introduction and Research Objectives In this chapter, the analysis methods that were used in this study and the results of this analysis are presented and discussed. The analyses are based on the objectives and research plan presented in Chapter 1. For convenience, the specific objectives of the study were included below as well. • Review the advantages and shortcomings of some of the existing frost and thaw depth models. • Develop accurate and reliable models for predicting the frost depth during freezing period in Michigan. • Develop a model to predict heave and the resulting pressure under the pavement or behind existing retaining structures due to freezing of frost-susceptible soils in Michigan. • Develop accurate and reliable model for predicting the thaw depth under the pavement in spring season in Michigan. • Investigate changes in pavement bearing capacity in the cycles of freeze and thaw; and • Develop a model to estimate when to begin and end SLR. To accomplish these objectives a comprehensive research plan was drawn as presented in Chapter 1 and depicted in Figure 4.1. Hypotheses Previous studies indicate that frost and thaw depths are a function of many variables including intensity and duration of the freezing period, water availability, soil permeability and capillarity, grain size and grain size distribution, and the soil thermal 85 conductivity. Hence, it was hypothesized that these variables are a function of the soil type such as clayey and sandy soils. Furthermore, it was hypothesized that the various missing soil parameters (such as insitu density, water content, grain size, grain size distribution, soil permeability and capillarity) can be expressed by one related property; the saturated thermal conductivity of the soil. Task 1 Conduct comprehensive literature review which includes the state of the art of modeling of freeze and thaw in soils and their applicability to this study, the state of art of modeling of frost heave and the resulting pressure behind retaining walls and the state of the practice of SHAs for forecasting frost and thaw depth and the time at which to post and remove SLR . Task 2 Develop heat transfer predictive model by firstly, scrutinize the existing models for frost, and thaw depth prediction using the MDOT database. Secondly, develop new statistical models for frost and thaw depth predictions using MDOT database. Thirdly, verify the accuracy of the developed model using MnDOT database. Task 3 Develop coupled heat and mass transfer models for prediction of frost heave and frost pressure by firstly estimate the freezing depth. Secondly, estimate the rate of flow of water to the frozen depth from a water supply. Thirdly, estimate the amount of heave pressure based on the frost heave model results. Finally, evaluate the accuracy of frost heave model using MDOT heave data. Task 4 Develop SLR policy based on different research on resilient modulus in literature and by using the results of thaw depth model Figure 4.1 Flow chart of the research plan 86 Frost Depth Existing frost depth prediction models can be classified into numerical, analytical, semiempirical, and empirical models. Some models require as inputs various thermal and hydraulic properties of soil and different meteorological data. Others require only the freezing index (FI) or the cumulative freezing degree day (CFDD). Any of these models can be used depending on the availability of the input data and the required accuracy. As mentioned in chapter 2, different numerical models were developed in the past three decades. In this chapter, UNSAT-H model was used to predict the frost depth. The results were then evaluated with the field data in Michigan. Furthermore, the accuracy of different analytical and semi-empirical frost depth prediction models including Stefan model (Jiji, 2009), Modified Berggren’s model (Aldrich et al., 1953) and Chisholm and Phang empirical model (Chisholm and Phang, 1983) were evaluated using the RWIS soil temperature data measured in the state of Michigan. Since none of the models yielded accurate results, revised empirical models that require only cumulative freezing degree day as input were developed. First, the data in the State of Michigan was used to develop an empirical model regardless of the soil types. Further, for model validation the 2003 to 2012 frost depth data from 8 stations in the State of Minnesota were used. Third, by considering the soil types two empirical models were developed for clayey and sandy soils. The two models were also evaluated using the Minnesota frost depth data. Finally, using the thermal conductivity data of each soil type, the two models were combined and one general model was developed which required CFDD and soil thermal conductivity. The accuracy of the general model was also checked using the frost depth data measured in the state of Minnesota. 87 4.3.1 UNSAT-H Model As stated in chapter 2, UNSAT-H is a one dimensional finite difference heat and mass balance model. The model inputs are meteorological data including air temperature, precipitation, solar radiation, wind speed, cloud cover, dew point, soil hydraulic, and thermal properties data. The soil can be modeled in different layers in UNSAT-H. For evaluating the UNSAT-H model, two of the RWIS sites located in the Lower Peninsula of Michigan were chosen. The pavement profile corresponding to both sites is illustrated in Figure 4.2. Hot mix asphalt (12.7 cm.) Gravel (30.5 cm.) Compacted Sand (61 cm.) Loose Sand (76 cm.) Figure 4.2 Schematic of the cross section of the two modeled pavement sites Table 4.1 shows the hydraulic and thermal properties that were used for each layer. Using these properties and the weather data as the boundary conditions, soil temperature profile was estimated. Figure 4.3 and Figure 4.4 show the estimated and the measured freezing depths as a function of time. It can be seen from the figures that the model predicted the freezing front up to three weeks earlier than the measured front. Although timing is not an issue in freezing period, it certainly is problematic during the thaw season. Also, in Cadillac station the model over predicted the frost depth by approximately 50 centimeters. One possible reason for the unfavorable results could be the fact that UNSAT-H does not consider latent heat of fusion in the 88 heat balance analysis. Therefore, while the model can predict temperatures above the freezing temperature quite accurately, the accuracy of temperature predictions below the freezing temperature is questionable. In addition, the use of the estimated input properties (such as thermal conductivity, volumetric specific heat, and hydraulic conductivity) contributed to the discrepancy between the estimated and the measured data. Table 4.1 Hydraulic and thermal properties for different layers in UNSAT-H model Asphalt 0.070 0.360 0.0050 5.60E-08 1.090 0.08257 3.9 Volumetric Specific Heat (kJ/(m3.oC) 2 Gravel 0.005 0.420 1.0000 1.00E+01 2.190 0.54338 1.25 1.36 0.020 0.375 0.0431 3.100 0.67742 1.5 2.39 Material Frost depth (cm) Sand θr 200 180 160 140 120 100 80 60 40 20 0 θs α Ks (cm/sec) n 4.63E-01 Thermal m (1-1/n) Conductivity (W/(m.oC)) Calculated frost depth (cm) Measured frost depth (cm) Date Figure 4.3 Calculated frost depth using UNSAT-H model and measured frost depth versus time in Waters station, Lower Peninsula, Michigan 89 Calculated frost depth (cm) Measured frost depth (cm) 200 180 Frost depth (cm) 160 140 120 100 80 60 40 20 0 Date Figure 4.4 Calculated frost depth using UNSAT-H model and measured frost depth versus time in Cadillac station, Lower Peninsula, Michigan 4.3.2 Freezing Index and Freezing Degree Day Calculation One of the common inputs to most analytical and semi-empirical models is the freezing index or the cumulative freezing degree day. Two different methods have been considered for calculating the cumulative freezing degree day; the Minnesota method (MnDOT, 2009) and Boyd method (Boyd, 1976). 4.3.2.1 Minnesota Cumulative Freezing Degree Day The cumulative freezing degree-day (CFDD) was calculated using Equation 4.1 (MnDOT, 2009): 𝑛 𝐶𝐹𝐷𝐷𝑛 = ∑ 𝐹𝑟𝑒𝑒𝑧𝑖𝑛𝑔 𝐷𝑒𝑔𝑟𝑒𝑒 𝐷𝑎𝑦 ≤ 0 𝑖=1 90 Equation 4.1 Freezing Degree day = ( Tmax + Tmin − 0o C) 2 Equation 4.2 Where 𝑇𝑚𝑎𝑥 = Maximum daily air temperature (°C); and 𝑇𝑚𝑖𝑛 = Minimum daily air temperature (°C). It should be noted that in the Minnesota method, the cumulative freezing degree day is reset on July 1 of each year and the freezing index is the maximum CFDD at the end of the winter season. Table 4.2 lists an example calculation of the cumulative freezing degree day. Table 4.2 Cumulative freezing degree day calculation, Waters station, Lower Peninsula Date Freezing Average Air Degree day Temperature (FDD) o o ( C) ( C-day) Cumulative freezing degree day (CFDD) (oC-day) Absolute Cumulative freezing degree day (CFDD) (oC-day) 11/16/2010 5.8 5.8 0 0 11/17/2010 3.2 3.2 0 0 11/18/2010 -0.95 -0.95 -0.95 0.95 11/19/2010 -0.75 -0.75 -1.7 1.7 11/20/2010 -1.25 -1.25 -2.95 2.95 11/21/2010 1.75 1.75 -1.2 1.2 11/22/2010 10.55 10.55 0 0 11/23/2010 4.15 4.15 0 0 11/24/2010 -2.3 -2.3 -2.3 2.3 11/25/2010 -2.35 -2.35 -4.65 4.65 11/26/2010 -6.4 -6.4 -11.05 11.05 11/27/2010 -1.95 -1.95 -13 13 11/28/2010 -1.15 -1.15 -14.15 14.15 11/29/2010 1.15 1.15 -13 13 4.3.2.2 Boyd Cumulative Freezing Degree Day If the CFDD is calculated and plotted as a function of time as shown in Figure 4.5, the graph will have a minimum value in the fall and a maximum value in spring. The Freezing Index 91 (FI) for that winter is estimated as the difference between the maximum and minimum cumulative degree days as shown in the Figure 4.5 (Boyd, 1976). Any spring or fall month that includes a seasonal maximum or a seasonal minimum degree days is called a “changeover” month. Boyd (1976) proposed Equation 4.3 that can be used for calculating the cumulative freezing degree days in the change-over month: 𝑌2 + 𝑁 ∗ 𝑋 ∗ 𝑌 = 𝑁2 𝑘2 Equation 4.3 where k = 2.5 constant; N= number of days in the month; 𝑋 = (𝑇 − 0℃ ); T= the average temperature in the change-over month of N days; and Y= Cumulative degree day of the change-over month. 1200 1000 CFDD (OC) 800 Freezing Index 600 400 200 0 Figure 4.5 Calculation of freezing index using cumulative freezing degree day 92 The solution of this equation yields two values for Y; a positive and a negative value. For the changeover month, the cumulative freezing degree days (CFDD) for the month can be calculated using Equation 4.4. Whereas, for all other months, the CFDD is calculated using Equation 4.5. 𝐶𝐹𝐷𝐷 = |𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒 𝑟𝑜𝑜𝑡 𝑜𝑓 𝑌| 𝐶𝐹𝐷𝐷 = 𝑁𝑋 Equation 4.4 Equation 4.5 Where all parameters are the same as before. CFDD values were calculated using both the Minnesota (Equation 4.1) and the Boyd methods and the data from different RWIS stations in Michigan. The results are plotted in Figure 4.6. It can be seen from the figure that the CFDD values calculated using Boyd equation differ from the CFDD values calculated using the Minnesota equation by more than 20 percent. This difference could be attributed to various reasons including: 1. The k value in Equation 4.3 is an empirical value based on the 10 years data collected at 22 stations across Canada. This value could change from one year to the next and from one region to another. 2. By using the monthly average temperature in the Boyd equation instead of the daily average temperature in calculating CFDD, the daily variations in the degree days are disregarded, which may lead to errors. Because of the above reasons, in this study, the CFDD values were calculated using the Minnesota equation (Equation 4.1). 93 Percent differnce with Boyd CFDD % 160 140 120 100 80 60 40 20 0 0 200 400 600 800 Minnesota CFDD (oC) 1000 1200 Figure 4.6 Comparison of calculated CFDD using Boyd (Boyd 1976) and Minnesota (MnDOT 2009) methods 4.3.3 Existing Frost Depth Prediction Models There are several frost prediction models that were developed. Some of these models are empirical in nature, some others are semi-empirical, and still others are mechanistically based. Some of these models are enumerated and discussed below. 4.3.3.1 Stefan’s Equation As stated before, Stefan solved the heat transfer phase-change problem for the special case of no heat transfer in the unfrozen zone (Jiji, 2009) and estimated the frost depths. His solution is one of the first frost depth prediction models (see Equation 2-19 of Chapter 2) and modified versions of his solution are still being used by some State Highway Agencies (SHAs). 94 In order to evaluate the accuracy of Stefan Equation relative to the measured data in Michigan, measured frost depths at different RWIS stations in the state of Michigan were used. Unfortunately the in-situ water content and dry density data of the soils were not available. Therefore, the soil water content and dry density were estimated using the graphs developed by the U.S Army Corps of Engineers (USACE, 1998) and shown in Figure 4.7. The various curves in the figure relate thermal conductivity to dry density and moisture content in the frozen and unfrozen conditions of various soil types. For each soil type, the measured thermal properties were used and its dry density and water content were estimated from the graphs. Next, the volumetric latent heat of fusion (L) was calculated using Equation 2-20; as it was expected, the calculated values of L decreased as the water content decreased. The values of the freezing index for the years 2010 and 2011 were calculated using the NOAA data obtained from the appropriate weather stations. Finally, Equation 2-19 was used to calculate the frost depths as a function of time for the two years. The details of frost depth calculation for each RWIS station were presented in appendix A. Figure 4.8 and Table 4.1 depicts the maximum calculated versus the maximum measured frost depth data for saturated condition. In Figure 4.8 the straight line is the line of equality between the measured and the calculated frost depth data. It can be seen from the figure that, for all soil types, the calculated maximum frost depths in saturated condition are much higher (more than 63.5 centimeters) than the measured values. The discrepancy between the measured and calculated data could be related to: 1. The volumetric heat capacity of the soil and water, which were not considered in Stefan’s Equation (Equation 2-19). 2. Errors in estimating the in-situ water content, dry density using the soil thermal conductivity and the Corps of Engineers graphs. 95 20 20 Dry Unit Weight (kN/m3) 19 17 16 1 2 .3 3 .4 4.5 5.7 6.8 8 9.1 19 7.9 Dry Unit Weight (kN/m3) K=(W/(m. C)) 14 K=(W/(m.oC)) 1.1 2.3 3.4 4.5 5.7 6.8 8 9.1 10 11 12 13 o 6.8 17 5.7 16 5.1 1.1 1.7 2.3 2.8 3.4 4 4.5 5.1 5.7 6.2 6.8 7.4 7.9 8.5 9.1 4.5 14 9.1 10 4.0 2.3 3.4 13 13 0 10 20 30 Water content (%) 0 40 10 20 30 Water content (%) Figure 4.7(a) Frozen sand and Gravel Figure 4.7(b) Unfrozen sand and Gravel 20 20 1.1 o K=(W/(m. C)) K=(W/(m.oC)) 2.3 2.3 19 2.8 Dry Unit Weight (kN/m3) 4.5 5.7 6.8 17 16 14 1 2 .3 3 .4 4.5 5.7 1.1 1.7 3.4 19 Dry Unit Weight (kN/m3) 40 3.4 4 4.5 17 5.1 5.7 16 14 6.8 8 13 13 0 10 20 30 Water content (%) 40 Figure 4.7(c) Frozen silt and clay 1.1 0 10 20 30 Water content (%) 40 Figure 4.7(d) Unfrozen silt and clay Figure 4.7 soil Thermal conductivity of different types of soil based on water content and dry density obtained by US army cold region and engineering laboratory (CRREL) (Edgar, 2014) 96 Line of equality between the measured and calculated data Maximum calculated frost depth (cm) 400 350 300 250 200 150 100 50 0 0 50 100 150 200 250 300 Maximum measured frost depth (cm) 350 400 Figure 4.8 The maximum frost depths predicted by Stefan equation versus the measured maximum frost depths in Michigan Given the significant differences between the measured frost depth data and the calculated ones using Equation 2-19, Stefan’s equation was abandoned and the modified Berggren’s equation was studied. The results are presented and discussed below. 4.3.3.2 Modified Berggren’s Equation Aldrich et al., (1953) made the two assumptions listed below and modified the original Berggren’s equation. Equation 2-21 was the resulting equation. 1. Heat transfer is one-dimensional problem and the soil is at its mean annual temperature before freezing begins (USACE, 1998). 97 Table 4.3 Maximum frost depth predicted by Stefan’s equation for RWIS stations Location Station Name Benzonia Cadillac Grayling Houghton Lake Lower Peninsula Ludington Reed City Waters Williamsburg Au Train Type of soil Year Loose Sand Dense Sand Dense Sand 2010-2011 2010-2011 2010-2011 Maximum Measured Frost Depth (cm) 86.4 116.8 157.5 Dense Sand 2010-2011 132.1 254.0 2010-2011 86.4 167.6 2010-2011 101.6 236.2 2010-2011 172.7 312.4 2010-2011 101.6 182.9 2009-2010 132.1 264.2 2010-2011 132.1 363.2 2009-2010 2010-2011 116.8 147.3 231.1 264.2 2009-2010 147.3 248.9 2010-2011 172.7 279.4 2010-2011 132.1 309.9 2010-2011 2010-2011 2010-2011 2010-2011 2010-2011 2010-2011 132.1 116.8 132.1 116.8 116.8 116.8 251.5 271.8 264.2 180.3 205.7 193.0 Loose Sand with clay Compacted Sand Loose Sand Compacted Sand Loose Sand Dense Sand Silty Clay Sand with Gravel and Silt Loose Sand Upper Peninsula Brevort Loose Sand Harvey Sand with Gravel and Silt Golden Lake Seney Cooks Michigamme St.Ignace Twin Lakes Engadine Dense Sand Dense Sand with Gravel Loose Sand Clayey Sand Clayey Sand Silty Clayey Sand Silty Clayey Sand Silty Clayey Sand 98 Maximum Calculated Stefan Eq. (cm) 221.0 271.8 256.5 2. At the beginning of the freezing season, the surface temperature decreases in a step-function manner from the mean annual temperature to some degrees below the freezing point and remains at this temperature (steady state) during the entire freezing season (Bianchini et al., 2012). In this study, the maximum frost depths were calculated using the modified Berggren’s equation for multilayered system (Equation 2-22).Table 4.4 shows an example of the step by step frost depth calculation using Equation 2-22. In the calculations, in order to obtain the correction factor (λ) from Figure 2.3, two dimensionless parameters (the thermal ratio (α) and the fusion parameter (μ)) must be calculated. As stated in Chapter 2, fusion parameter (μ) depends on the volumetric heat capacity and latent heat of fusion of every layer and was calculated using Equation 4.6. The thermal ratio (α) is a fixed number for all layers and depends on the FI and annual average temperature; it was calculated using Equation 4.7. 𝜇= 𝐶 ∗𝑣 𝐿 𝑠 𝛼= 𝑣0 𝑣𝑠 Equation 4.6 Equation 4.7 Where C=volumetric heat capacity (KJ/m3); L= latent heat of fusion (KJ/m3); vs= average temperature differential= n(FI)/t; t= duration of winter period (used in calculation of vs); v0= initial temperature differential = annual average temperature -0 oC; and All other parameters are as before. After calculating α and μ, the correction factor λ was obtained from Figure 2-3 and the maximum frost depth was calculated using Equation 2-22 (See Table 4.4). 99 Using the Modified Berggren’s equation, the maximum frost depths were calculated for the RWIS stations in Michigan. The details of frost depth calculation for each station were presented in appendix A. The maximum calculated frost depths were compared to the maximum measured frost depth data in the State of Michigan in Figure 4.9 and Table 4.5. It can be seen from the figure that the Modified Berggren’s equation leads to more accurate results than the Stefan equation. However, the differences between the calculated and measured values in some cases are more than 51 centimeters. The discrepancy between the measured and calculated data could be related to: 1. Equation 2-22 does not account for the water movement in the soil. 2. Potential errors in estimating the thermal conductivity, water content and dry density. Given the substantial differences between the measured frost depth data and the calculated ones using Equation 2-22, the Modified Berggren’s equation was also abandoned and the Chisholm and Phang equation was studied. The results are presented and discussed in the next section. 4.3.3.3 Chisholm’ and Phang’s Equation One of the first empirical equations, which relate CFDD and frost depth, was developed by Chisholm and Phang in 1980 to predict frost depths under asphalt pavements in Ontario, Canada (Equation 2-23). It should be noted that since the local daily air temperature data were not available on a daily basis, the CFDD values were calculated using Boyd approach (Chisholm and Phang, 1980). Their results indicated that the frost depth predictions were within 30.5 centimeters of the measured values. 100 Table 4.4 Frost depth calculation using the modified Berggren’s equation, Benzonia, Lower Peninsula, Michigan Column No. 1 2 3 4 5 Material γd w d C k HMA 21.7 20.4 0 0.13 894 0.075 0.30 1118 1.49 2.60 6 7 Ĺ= ∑Ld/∑d 14792 0 50300 35508 L 8 9 Ĉ= ∑Cd/∑d 894 1043 μ= Ĉ / Ĺ *vs 0 10 11 12 13 14 15 λ R ∑R ∑R+R/2 FI ∑FI 0 0.3 0.4 0.0 0.3 0.1 0.5 0 82 0 82 Gravel 0.17 0.56 Loose 0.16 0.58 0.8 0.7 19.6 0.09 0.58 1565 2.42 57938 48362 1341 1.0 372 457 Sand γd = unit weight (kN/m3); w = water content (%); d = layer depth (m); C= volumetric heat capacity(KJ/m3); k= thermal conductivity (W/( m.oC ); L= latent heat of fusion (KJ/m3); μ= fusion parameter; λ= correction factor; R= thermal diffusivity ((m2.oC)/W); FI=freezing Index (oC-day) Step by step calculation 1. k values were obtained from the laboratory measurements (Table 3.4). γd and w values were obtained from Figure 4.7 using k values. d values were obtained from the pavement profile of RWIS station, and C values are assumed based on the soil type (Columns 1-5). 2. L and μ were calculated using Equation 2.20 and 4.6, respectively (Columns 6 and 9). 3. α can be calculated by Equation 4.7 as follow: FI=445 (oC-day) ; vs=0.9(445)/t = 3.1 ; v0= 9.2-0=9.2  α=v0/vs=2.97 4. Using α and μ values, λ can be obtained from Figure 2.3 for each layer (Columns 10) 5. R values were calculated as R= d/K for each layer (Columns 11) 6. Freezing index required for each layer to freeze were calculated using Equation 2.22 (Columns 14) 7. The summation of freezing indexes in column 14 should be approximately equal to the seasonal freezing index (FI=445) 101 Line of equality between the measured and calculated data Maximum calculated frost depth (cm) 400 350 300 250 200 150 100 50 0 0 50 100 150 200 250 300 Maximum measured frost depth (cm) 350 400 Figure 4.9 Measured maximum frost depths in Michigan versus the maximum calculated ones using the modified Berggren’s equation Equation 2-23 was used to predict the maximum monthly measured frost depths at different RWIS stations in the state of Michigan. Figure 4.10 shows the results. It can be seen that in most cases, Equation 2-23 underestimates the maximum monthly frost depths. In fact, in some cases for small values of CFDD, the calculated frost depths could be negative. The differences between the predicted and the measured frost depths could be as high as 76 centimeters. Based on the fact that the empirical equation was developed using the measured frost depth data in Ontario, it can be concluded that the equation is regional and calibration is required for using it in other regions. 102 Table 4.5 Maximum frost depth predicted by Modified Berggren’s equation for RWIS stations Location Station Name Benzonia Cadillac Grayling Houghton Lake Lower Peninsula Ludington Reed City Waters Williamsburg Au Train Type of soil Year Loose Sand Dense Sand Dense Sand 2010-2011 2010-2011 2010-2011 Maximum Measured Frost Depth (cm) 86.4 116.8 157.5 Dense Sand 2010-2011 132.1 2010-2011 86.4 2010-2011 101.6 2010-2011 172.7 2010-2011 101.6 2009-2010 132.1 2010-2011 132.1 2009-2010 2010-2011 116.8 147.3 2009-2010 147.3 2010-2011 172.7 2010-2011 132.1 2010-2011 2010-2011 2010-2011 2010-2011 2010-2011 2010-2011 132.1 116.8 132.1 116.8 116.8 116.8 Loose Sand with clay Compacted Sand Loose Sand Compacted Sand Loose Sand Dense Sand Silty Clay Sand with Gravel and Silt Loose Sand Upper Peninsula Brevort Loose Sand Harvey Sand with Gravel and Silt Golden Lake Seney Cooks Michigamme St.Ignace Twin Lakes Engadine Dense Sand Dense Sand with Gravel Loose Sand Clayey Sand Clayey Sand Silty Clayey Sand Silty Clayey Sand Silty Clayey Sand 103 Maximum Calculated Stefan Eq. (cm) 101.6 152.4 162.6 160.0 96.5 139.7 170.2 109.2 149.9 213.4 154.9 175.3 132.1 154.9 215.9 162.6 177.8 193.0 129.5 148.9 132.1 Line of equality between the measured and calculated data Maximum Calculated Frost depth (cm) 400 350 300 250 200 150 100 50 0 0 50 100 150 200 250 300 350 400 -50 Maximum Measured Frost Depth (cm) Figure 4.10 Measured maximum frost depths versus calculated ones using Chisholm’ and Phang’s equation 4.3.4 Frost Depth Empirical Models Since none of the existing models yielded accurate frost depth results, new empirical models were developed in this study using the RWIS data in the State of Michigan. These new models are presented below. Among the 25 RWIS stations located in the Upper and Lower Peninsulas of Michigan, the data from 18 stations were used for developing empirical models. Seven stations were not considered due to incomplete data (some of the data are missing). For air temperature, the databases from the nearest NOAA stations were obtained and used due to higher consistency and accuracy relative to the RWIS database. Sites not situated close to NOAA stations are not considered in this study. Due to data availability, only data from the winter of 2011 was used, 104 except for Au Train, Brevort and Harvey sites where data from the 2010 winter were also available. First, the air temperature data were used to calculate the CFDD for each RWIS station location. Second, the measured frost depth data in all RWIS stations and the calculated CFDD values were used to develop one simple statistical prediction equation regardless of the soil type. This resulted in Equation 4.8. 𝑃 = 4.759 (𝐶𝐹𝐷𝐷)0.5339 Equation 4.8 Equation 4.8 is more or less similar to Equation 4.9 , which was developed by the U.S. Corps of Engineers (Yoder, 1975) 𝑃 = 4.210 (1.8 𝐶𝐹𝐷𝐷)0.478 Equation 4.9 Where P = frost depth (cm); and CFDD = cumulative freezing degree day (oC – day). The results of both equations are depicted in Figure 4.11. The data in the figure indicate that Equation 4.9 under predicts the majority of the data. On the other hand, Equation 4.8 represents the measured frost depths more accurately. In fact, the calculated coefficient of determination (R2) is 0.91. Figure 4.12 show the calculated frost depths using Equation 4.8 versus the measured ones. The solid straight line is the locus of equality between the measured and calculated data. It can be seen that, the majority of the calculated frost depth data are within a few centimeters from the measured values and the maximum difference is 25-centimeter. As stated in Chapter 3, frost depth data measured from 2003 to 2012 in 8 stations in the State of Minnesota were requested and received from MnDOT. Equation 4.8 was then used to calculate the frost depth data at all 8 stations and for the ten year period. The measured frost depth data and the calculated ones are depicted in Figure 4.13. The straight line in the figure 105 indicates the line of equality between the measured and the calculated values. Examination of the data shown in the figure indicates that Equation 4.8 does not predict the measured frost depth data in Minnesota accurately. In fact, Equation 4.8 over predicts the Minnesota frost depth data by as much as 102 centimeters. Moreover, the calculated coefficient of determination (R2) for the Minnesota data is 0.77, which is much lower than the calculated R2 of 0.91 for the Michigan data. Other performance metrics of Equation 4.8 and Equation 4.9 are shown in Table 4.6. In this table, standard error of the estimate (SEE), mean absolute error (MAE) and mean absolute percentage error (MAPE) are shown for each model. The SEE is proportional to the width of the confidence interval so smaller SEE is an indication of a better fit. As can be seen in, Table 4.6 for the Michigan, the SEE value of Equation 4.9 is 60% larger than SEE of Equation 4.8. On the other hand, Equation 4.9 performs better in the State of Minnesota. In general, all of the performance metrics indicate that Equation 4.8 does not predict frost depth data in the State of Minnesota accurately. U.S.Corp of Engineers 200 Clayey soil P = 4.7591CFDD0.5339 R² = 0.9135 n= 158 160 Frost depth (cm) Sandy soil 120 80 40 0 0 200 400 Cumulative freezing degree days (oC-day) 600 Figure 4.11 Frost depths versus cumulative freezing degree day for clayey and sandy soils in the State of Michigan 106 Line of equality between the measured and calculated data Clayey soil-Michigan Sandy soil-Michigan Calculated frost depth (cm) 200 n= 158 160 120 80 40 0 0 40 80 120 Measured frost depth (cm) 160 200 Figure 4.12 Measured frost depths in State of Michigan versus the calculated ones using Equation 4.8 Calculated frost depth (cm) Line of equality between the measured and calculated data Calyey soil-Minnesota Sandy soil-Minnesota 240 n= 621 200 160 120 80 40 0 0 40 80 120 160 Measured frost depth (cm) 200 240 Figure 4.13 Measured frost depth in State of Minnesota versus the calculated ones using Equation 4.8 107 Table 4.6 Performance metrics of Equation 4.8 and Equation 4.9 for frost depth estimations in the State of Michigan and Minnesota Location Michigan Minnesota Model Standard Error of the Estimate (SEE) Mean Absolute Error (MAE) Equation 4.8 15.0 11.2 Mean Absolute Percentage Error (MAPE) 0.18 Equation 4.9 19.1 12.7 0.18 Equation 4.8 24.4 18.3 0.29 Equation 4.9 17.5 13.2 0.22 Figure 4.14 and Figure 4.15 show the calculated frost depths using Equation 4.9 versus the measured ones in the state of Michigan and Minnesota, respectively. The solid straight line in the figures is the locus of equality between the measured and the calculated data. It can be seen that Equation 4.9 predicts frost depths in clayey soils better than in sandy soils in both states. For the latter soils, the differences between the measured and calculated values could be as high as 63.5-centimeter. Please note that, like Equation 4.8, Equation 4.9 does not separate sandy from clayey soils. The results shown in Figure 4.11 were further scrutinized to improve the accuracy of Equation 4.8. Previous studies indicate that frost depths are a function of many variables including intensity and duration of the freezing period, water availability, soil permeability and capillarity, grain size and grain size distribution, and the soil thermal conductivity. Hence, it was hypothesized that these variables are a function of the soil type such as clayey and sandy soils. Therefore, the frost depth data measured at various RWIS stations in the state of Michigan was divided into two groups according to soil type at the stations; clayey and sandy soils. It should be noted that dividing the data into two groups of clayey and sandy soil was based on the soil log provided by MDOT for each station. 108 Line of equality between the measured and calculated data Clayey soil-Michigan Sandy soil-Michigan Calculated frost depth (cm) 200 160 n=158 120 80 40 0 0 40 80 120 Measured frost depth (cm) 160 200 Figure 4.14 Measured frost depth in State of Michigan versus the calculated ones using Equation 4.9 Calculated frost depth (cm) Line of equality between the measured and calculated data Calyey soil-Minnesota Sandy soil-Minnesota 240 200 n= 621 160 120 80 40 0 0 40 80 120 160 Measured frost depth (cm) 200 240 Figure 4.15 Measured frost depth in State of Minnesota versus the calculated ones using Equation 4.9 109 After dividing the data per soil type, a mathematical power function was used to model each group of frost depth data as a function of the calculated CFDD. This resulted in Equation 4.10 Equation 4.11 for clayey and sandy soils, respectively. For Clayey Soils 𝑃 = 5.3858 (𝐶𝐹𝐷𝐷)0.4896 Equation 4.10 For Sandy Soils 𝑃 = 4.6473 (𝐶𝐹𝐷𝐷)0.5423 Equation 4.11 Figure 4.16 and Figure 4.17 depict the measured frost depth data in clayey soils in Michigan and the calculate frost depth data using Equation 4.10. Figure 4.16 also show the U.S. Corps of Engineers equation (Equation 4.9). It can be seen from the figure that Equation 4.10 and Equation 4.9 fit the data very well. The coefficient of determination for both equations is 0.94. Further, Figure 4.17 depicts the measured frost depth data in clayey soils in Michigan versus the frost depth data calculated using Equation 4.10. The solid line in the figure is the line of equality between the measured and the calculated data. The results in the figure indicate that Equation 4.10 predict the frost depth data in clayey soils in Michigan more accurately. Figure 4.18 Figure 4.19 depict the measured frost depth data in sandy soils in Michigan and the calculate frost depth data using Equation 4.11. Figure 4.18 also show the U.S. Corps of Engineers equation (Equation 4.9). It can be seen from the figure that Equation 4.11 represents the measured frost depth data much better than Equation 4.9. Indeed, the coefficient of determination is 0.91 and 0.76 for Equation 4.11 and Equation 4.9, respectively. Further, Figure 4.19 depicts the measured frost depth data in clayey soils in Michigan versus the calculated values using Equation 4.11. The solid line in the figure is the line of equality between the measured and the calculated data. The results in the figure indicate that Equation 4.11 provides similar prediction of the measured frost depth data to Equation 4.8. Since 110 the majority of the data were from sandy soil, the similarity between the predictions of both equations was expected. U.S.Corp of Engineers Frost depth (cm) 200 P = 5.3858(CFDD)0.4896 R² = 0.94 n=29 160 120 80 40 0 0 200 400 Cumulative freezing degree days (oC- day) 600 Figure 4.16 Measured frost depths in Michigan versus cumulative freezing degree day for clayey soil showing the best fit and the U.S. Corps of Engineers equations Calculated frost depth (cm) Line of equality between measured and calculated data Clayey soil-Michigan 200 160 n= 29 120 80 40 0 0 40 80 120 160 Measured frost depth (cm) 200 Figure 4.17 Measured versus calculated frost depth data in clayey soils in Michigan (Equation 4.10) 111 U.S.Corp of Engineers 200 P = 4.6473(CFDD)0.5423 R² = 0.9119 n= 129 Frost depth (cm) 160 120 80 40 0 0 200 400 600 Cumulative freezing degree days (oC-day) Figure 4.18 Frost depths versus cumulative freezing degree day for sandy soil showing the best fit and the U.S. Corps of Engineers equations in the State of Michigan Line of equality between the measured and calculated data Calculated frost depth (cm) Sandy soil-Michigan 200 160 n= 129 120 80 40 0 0 40 80 120 160 Measured frost depth (cm) 200 Figure 4.19 Measured versus calculated frost depths in sandy soils in Michigan (Equation 4.11) It should be noted that for frost depth more than 127-centimeter in sandy soils in both states, Equation 4.11 underestimates the measured data by as much as 25-centimeter. 112 Examination of the results depicted in Figure 4.16 through Figure 4.19 indicates that Equation 4.10 predicts the frost depth data in clayey soils better than Equation 4.11 in sandy soils. The main reason is that the variability of the measured frost depth data in sandy soils is higher than that in clayey soil. Such variability is a function of the grain size and grain size distribution, which impact the distribution of water and the hydraulic conductivity of the soils. Unfortunately, such data are not available at this time to improve Equation 4.11. Nevertheless, the equation does predict the frost depth data in sandy soils relatively accurately. One important point should be noted is the number of measured data points in clayey soils is much less than in sandy soils. A total of 29 data points are available for clayey soils, whereas 129 data points are available in sandy soils. Once again, to evaluate the accuracy of Equation 4.10 and Equation 4.11 , they were used to predict the measured frost depths in clayey and sandy soils in the State of Minnesota. Figure 4.20 and Figure 4.21 depict the results. It should be noted that for clayey soil (Figure 4.20) the number of measured data points is 374 while for sandy soil (Figure 4.21) it is 247. Examinations of Figure 4.20 and Figure 4.21 indicates that the prediction of frost depth data in clayey and sandy soils in Minnesota using Equation 4.10 and Equation 4.11 is better and more accurate than the prediction using Equation 4.8 (see Figure 4.13, Figure 4.20, and Figure 4.21). In fact, for clayey and sandy soil the calculated coefficients of determination (R2) are 0.88 and 0.9, respectively. The relatively high values of R2 indicate that Equation 4.10 and Equation 4.11 predict the frost depth data in clayey and sandy soils in the states of Michigan and Minnesota relatively accurately. In order to further evaluate the performance of Equation 4.8 to Equation 4.11, the performance metrics of all four equations in clayey and sandy soils for both states listed in Table 4.7 were examined. 113 Line of equality between the measured and calculated data Clayey soil- Minnesota Calculated frost depth (cm) 240 n= 374 200 160 120 80 40 0 0 40 80 120 160 Measured frost depth (cm) 200 240 Figure 4.20 Measured frost depths in clayey soil in the state of Minnesota versus the frost depth values calculated using Equation 4.10 Line of equality between the measured and calculated data Sandy soil- Minnesota Calculated frost depth (cm) 240 n= 247 200 160 120 80 40 0 0 40 80 120 160 Measured frost depth (cm) 200 240 Figure 4.21 Measured frost depths in sandy soil in the state of Minnesota versus the frost depth values calculated using Equation 4.10 114 As expected, the performance of Equation 4.9 and 4.10 are almost the same in both states. On the other hand, while performance of Equation 4.10 in Minnesota is not as good as in Michigan, it is substantially better than Equation 4.8 in both states. In addition, as expected, since the majority of the data in Michigan were from sandy soil, Equation 4.11 and Equation 4.8 performance are similar in both states. Table 4.7 Performance metrics of Equation 4.8 and Equation 4.9 for frost depth estimations in the State of Michigan and Minnesota Location Soil type Clayey Soil Michigan Sandy Soil Clayey Soil Minnesota Sandy Soil Model Standard Error of the Estimate (SEE) Mean Absolute Error (MAE) Equation 4.8 11.9 8.4 Mean Absolute Percentage Error (MAPE) 0.15 Equation 4.9 7.4 5.1 0.11 Equation 4.10 7.1 4.6 0.11 Equation 4.8 15.5 11.7 0.18 Equation 4.9 20.8 14.2 0.19 Equation 4.11 15.2 11.7 0.18 Equation 4.8 28.4 21.8 0.35 Equation 4.9 13.7 11.2 0.23 Equation 4.10 15.0 11.9 0.25 Equation 4.8 15.7 12.4 0.17 Equation 4.9 22.1 16.5 0.18 Equation 4.11 17.3 13.0 0.18 The thermal conductivity of any soil type (see Chapter 2) depends upon its water content, dry density, void distribution, and grain size and grain size distribution. These physical properties vary substantially from one soil type to another. Therefore, the disturbed clayey and sandy soil samples that were obtained by MDOT from various RWIS stations were saturated and the thermal conductivity of each soil type was measured in the laboratory at Michigan State 115 University using the KD2 pro. The results are listed in table 3-4 of Chapter 3. To consolidate Equation 4.10 and Equation 4.11 into one equation, it was hypothesized that the various missing soil parameters (such as insitu density, water content, grain size, grain size distribution, soil permeability and capillarity) can be expressed by one related property; the saturated thermal conductivity of the soil. Based on the hypothesis, the statistical parameters of Equation 4.10 and Equation 4.11 were correlated to the average thermal conductivity of each soil type. The statistical parameters of the two equations were then replaced by the resulting correlation equation, which yielded Equation 4.12 for both clayey and sandy soils. Figure 4.22 and Figure 4.23 show the correlation between the statistical parameters and the average thermal conductivity of each soil type. Unfortunately, only two data points (two soil types) were available, hence the best correlation between the statistical parameters and the average thermal conductivity is a straight line as shown in the figures. It should be noted that such straight line correlations may not be accurate and may result in errors in the resulting frost prediction equations. To produce more accurate nonlinear equations (power, exponential or logarithmic function), data from three or more soil types must be available. Unfortunately, this was not the case and the straight line equations are the best scenario for the given data. Nevertheless, the equations for the two straight lines in Figure 4.22 and Figure 4.23 were used to replace the statistical constants of Equation 4.10 and Equation 4.11. Equation 4.12 is the resulting equation for both types of soils clayey and sandy. 𝑃 = (−0.7392 𝑘 + 6.4408) ∗ (𝐶𝐹𝐷𝐷)(0.0035𝑘+0.4846) Where k= the average thermal conductivity of the soil (W/(m.oC)); and All other parameters are the same as before. 116 Equation 4.12 0.55 0.54 b = 0.0035k + 0.4846 0.53 b 0.52 0.51 0.5 0.49 0.48 0 0.5 1 1.5 2 2.5 o Thermal conductivity (W/(m. C)) 3 Figure 4.22 Correlation between the statistical power coefficient (b) of Equation 4.10 and a Equation 4.11 and the corresponding average thermal conductivity of the soil 5.5 5.4 5.3 5.2 5.1 5 4.9 4.8 4.7 4.6 a = -0.7392k + 6.4408 0 1 2 Thermal conductivity (W/(m. oC)) 3 Figure 4.23 Correlation between the statistical coefficient (a) of Equation 4.10 and Equation 4.11 and the corresponding average thermal conductivity of the soil Equation 4.12 was then used to calculate the frost depths in clayey soils in the States of Michigan and Minnesota. The inputs to the equation consisted of the calculated CFDD for each state and the average measured thermal conductivity of the soil samples obtained from MDOT. 117 Figure 4.24 and Figure 4.25 depict the calculated and the measured frost depths in Michigan and in Minnesota, respectively. Comparing the results shown in the two figures and those shown in Figure 4.17 and Figure 4.20 using Equation 4.10 indicate that the two equations produce similar results for clayey soils. Similarly, for sandy soils in Michigan and Minnesota, Equation 4.11 and Equation 4.12 produced almost the same results. These results can be seen in Figure 4.19 and Figure 4.21 for Equation 4.11 and in Figure 4.24 and Figure 4.50 for Equation 4.12. Performance metrics of Equation 4.12 in clayey and sandy soils for both states are listed in Table 4.8. As expected, all performance metrics indicate the similarity between the performance of Equation 4.10, 4.11, and 4.12. Line of equality between the measured and calculated data Clayey soil- Michigan Sandy soil-Michigan Calculated frost depth (cm) 200 n = 158 160 120 80 40 0 0 40 80 120 Measured frost depth (cm) 160 200 Figure 4.24 Calculated frost depths using Equation 4.12 versus the measured frost depth in clayey and sandy soil in the State of Michigan 118 Line of equality between the measured and calculated data Clayey soil-Minnesota Sandy soil-Minnesota Calculated frost depth (cm) 240 n= 621 200 160 120 80 40 0 0 40 80 120 160 Measured frost depth (cm) 200 240 Figure 4.25 Calculated frost depths using Equation 4.12 versus the measured frost depth in clayey and sandy soil in the State of Minnesota Table 4.8 Performance metrics of Equation 4.12 for frost depth estimations in the State of Michigan and Minnesota Location Soil type Model Michigan Clayey Soil Sandy Soil Equation 4.12 Minnesota Clayey Soil Sandy Soil Equation 4.12 Standard Error of the Estimate (SEE) Mean Absolute Error (MAE) 7.112 4.826 Mean Absolute Percentage Error (MAPE) 0.11 17.272 11.938 0.18 14.986 11.938 0.25 16.002 12.446 0.16 Recall that the average thermal conductivity values obtained from seven different soil samples (two soil types) were used in the prediction of the frost depths in Michigan and Minnesota. This implies that: 119 1. If thermal conductivity data of more soil types are available, the prediction of frost depth could improve. 2. Equation 4.12 could perhaps be used at the regional level to estimate the frost depth data. To further evaluate the validity of Equation 4.12, a statistical model was developed using the measured frost depth and the calculated CFDD for both soil types in Minnesota. The results are shown in Figure 4.26 and Figure 4.27, respectively. The dashed and solid curves in Figure 4.26 and Figure 4.27 represent the statistical model and Equation 4.12, respectively. Examination of Figure 4.26 indicates that Equation 4.12 and the statistical model produce almost the same results for clayey soils in Minnesota. On the other hand, the data in Figure 4.27 indicate that for sandy soil the differences between the results of the two models are less than 12.5 centimeters. This implies that Equation 4.12 could be used in Minnesota without calibration. Statistical power function Clayey soil-Minnesota Proposed model (Eq.4-12) 240 Frost depth (cm) 200 160 n= 374 120 80 40 0 0 400 800 1200 1600 o Cumulative freezing degree days ( C-day) Figure 4.26 Frost depths versus cumulative freezing degree-day for clayey soil showing the best fit statistical model and the proposed model (Equation 4.12) in Minnesota 120 Statistical power function Sandy soil-Minnesota Poposed model (Eq.4-12) 280 Frost Depth (cm) 240 200 n= 247 160 120 80 40 0 0 400 800 1200 o Cumulative freezing degree days ( C-day) Figure 4.27 Frost depths versus cumulative freezing degree-day for sandy soil showing the best fit and proposed model (Equation 4.12) in the State of Minnesota Frost Heave Frost heave refers to the uplifting of the ground surface caused by freezing of water within the soil layers. In cold regions, frost heave could cause uplifting of the pavement structure, shoulders, and even unprotected foundations of bridges and trusses supporting highway signs and utility lines (Liu et al., 2012). Frost heave can be influenced by various conditions including: 1. Frost Susceptibility – In general, the frost susceptibility of a soil is a function of its grain size and grain size distribution, which affect its capillarity and hydraulic conductivity (ACPA, 2008). There are various methods and criteria for the determination of soil frost susceptibility. In general, frost susceptibility could be affected by soil type. In coarse material such as gravel and coarse sand hydraulic conductivity is high and capillary potential is low, whereas clay has low hydraulic conductivity and high capillary potential. Only fine sand and 121 silt seem to have a balance between hydraulic conductivity and capillary potential. Figure 4.28 illustrates the dual effect of hydraulic conductivity and capillary potential on frost susceptibility. One of the most common criteria regarding frost susceptibility is based on the grain size distribution and the percent passing sieve number 200. Figure 4.29 and Table 4.9 show the susceptibility criteria developed by the U.S Corp of Engineers (COE). The Canadian Department of Transportation developed another soil frost susceptibility criterion that also based on soil grain sized distribution as shown in Figure 4.30. Figure 4.28 Effect of capillary and permeability on frost susceptibility (ACPA, 2008) 122 Figure 4.29 Heaving Rate in laboratory test on different disturbed soil types (COE, 1984) 123 Table 4.9 Frost susceptibility classification (COE, 1984) Typical Soil Types Under Unified Soil Classification System Frost Group Soil Percentage Finer Than 0.02 mm by Weight Non-frost susceptible Gravel Crushed stone Crushed rock 0- 1.5 GW, GP Sands 0- 3 SW,SP Gravel Crushed stone Crushed rock 1.5- 3 GW,GP 3- 10 SW,SP Possibly frost susceptible, requires lab tests Sands S1 Gravely soils 3- 6 S2 Sandy soils 3- 6 F1 Gravely soils 6- 10 Gravely soils 10- 20 Sands 6- 15 F2 F3 F4 GW, GP, GWGM, GP-GM SW, SP, SW-SM, SP-SM GM , GW-GM, GP-GM GM , GW-GM, GP-GM SM , SW-SM, SP-SM Gravely soils Over 20 GM, GC Sands, except very fine silty sands Over 15 SM, SC Clays, PI>12 --- CL, CH Silts --- ML, MH Very fine, silty sand Over 15 SM Clays, PI< 12 --- CL, CL-ML --- CL, ML and SM, CL, CH and ML, CL, CH, ML and SM Varved clays and other fine-grained, banded sediments 124 Figure 4.30 Frost susceptibility criteria, Canadian Department of Transportation (Edgar, 2014) 2. Below Freezing Temperature - As stated in Chapter 2, freezing point depression occurs in pore water because of different reason such as intermolecular forces between water and soil (soil water surface tension) and salt solution. Therefore, pore water starts to freeze when the air temperatures and consequently the ground surface temperature drops below the freezing temperature of 0oC. The rate of water freezing is a function of the actual temperature below freezing and its duration. Colder and more sustainable below freezing temperatures accelerate the freezing rate and increases the depth of frost penetration and consequently increases ground heave. Snow cover acts like insulator reducing frost depth substantially unless the air temperature and consequently the soil surface temperature drop significantly below the freezing temperature. However, for safety reasons, snow is typically removed from 125 the pavement surface and accumulated near the shoulder as soon as possible. This causes higher frost depth and higher frost heave under the pavements relative to other areas covered by snow (Yoder, 1975). Further, salt and other deicing chemicals (typically used on roads during winter season) decrease the temperature at which water starts to freeze and causes decreases in frost depth and frost heave. 3. Availability of Water Source – If no free water is available, no water frost action will take place, hence, a source of water should be available under the pavement to start the free water freezing process. The water source could be as deep as 6 meters (Edgar, 2014). If the ground water level is shallow, frost heave can be observed even in course material (COE, 1984). Figure 4.31 shows the ASSHTO four different environmental regions in the United States. Only two regions, wet-freeze, and dry freeze are subjected to water freezing under the pavements. The wet freeze region is considered to be the most frost susceptible region (ACPA, 2008). As can be seen, the state of Michigan is located in the most frost susceptible region, the wet-freeze region. Hence, the estimation of frost depths and frost heave are two important factors that are typically considered in the design of pavement and bridge and other structural foundations. 4. The Presence or Absence of Insulation- See chapter 2 for frost heave mitigation methods and the impacts of insulation on frost depth. 4.4.1 Insulation Effect on the Frost Depth As stated before, insulation is typically used to reduce heat flow and prevent heat loss in temperatures below 0 oC. Since the thermal conductivity of the insulation material is low, the heat loss decreases across the soil layer and the temperature remains above freezing point. Equation 4.13 that is based on conservation of energy governs the Temperature variation in a soil 126 layer (Jiji, 2009). Since no energy is generated in the freezing process, the equation can be rewritten as Equation 4.14. Using Fourier’s law, the heat flux can then be calculated using Equation 4.15. Figure 4.31 AASHTO four environmental regions (ACPA, 2008) Equation 4.13 𝐸̇𝑖𝑛 + 𝐸̇𝑔 − 𝐸̇𝑜𝑢𝑡 = 𝐸̇ 𝐸̇𝑖𝑛 − 𝐸̇𝑜𝑢𝑡 = 𝐸̇ 𝐸̇𝑖𝑛 − 𝐸̇𝑜𝑢𝑡 = −𝑘 Where 𝜕 2𝑇 𝜕𝑧 2 𝐸̇ = rate of energy change within the region (J/s); 𝐸̇𝑖𝑛 = rate of energy added (J/s); 𝐸̇𝑔 = rate of energy generated (J/s); 𝐸̇𝑜𝑢𝑡 = rate of energy removed (J/s); T= temperature in the soil layer (oC); k= Thermal conductivity of the soil layer (W/(m.oC)); and 127 Equation 4.14 Equation 4.15 z= depth from the ground surface (m). The rate of energy change within the region can be calculated using Equation 4.16. Thus, Equation 4.14 can be rewritten as Equation 4.17 (Edgar, 2014). Ė = 𝐶 ∂T ∂t 𝜕𝑇 𝜕 2𝑇 𝐶 = −𝑘 2 𝜕𝑡 𝜕𝑧 Equation 4.16 Equation 4.17 It should be noted that Equation 4.17 does not consider the phase change effect in the soil layer but since in RPF latent heat of fusion is negligible, this equation can be used for modeling an insulation layer. By assuming that surface temperature varies in a sinusoidal manner, solution to Equation 4.17 can be obtained using Equation 4.18 (Edgar, 2014). −𝑧 𝑧 𝑇(𝑧, 𝑡) = 𝑇𝑎𝑣𝑒 + 𝐴0 ∗ 𝑒 𝑑 ∗ sin (𝜔𝑡 − ) 𝑑 Equation 4.18 Where 𝑇(𝑧, 𝑡)= temperature variation at depth for each time interval (oC); Tave= the average temperature in soil layer (oC); A0= amplitude of the sine wave which relates to surface temperature fluctuation; d= depth that relates the reduction in temperature fluctuation A0 to depth (m); ω = time frequency; and All other parameters are the same as before. As stated above, the parameter “d” in Equation 4.18 is the characteristic depth that relates the reduction in surface temperature fluctuation to depth and can be calculated using Equation 4.19 (Edgar, 2014). 𝑑=( 2𝑘 1 )2 𝐶𝜔 Where all parameters are the same as before. 128 Equation 4.19 According to Equation 4.19, adding a low thermal conductivity layer could significantly influence the temperature pattern in the soil. In fact, adding an RPF layer has the same effect as adding additional soil to the layer. Therefore, the thickness of the RPF can be modeled as 𝑡𝑅𝑃𝐹 𝑘𝑅𝑃𝐹 = 𝑑 𝑠𝑜𝑖𝑙 ∗ √ 𝑘𝑠𝑜𝑖𝑙 Equation 4.20 Where tRPF = insulation thickness (cm); dRPF = depth of the soil layer (cm); 𝑘𝑠𝑜𝑖𝑙 = thermal conductivity of soil (W/(m.oC)); and 𝑘𝑅𝑃𝐹 = the effective thermal conductivity of insulation layer (W/(m.oC)). The thickness of the RPF can be calculated using Equation 4.21 𝑡𝑅𝑃𝐹 = (((−0.7392 𝑘 + 6.4408) ∗ (𝐶𝐹𝐷𝐷) (.0035𝑘+0.4846) 𝑘𝑅𝑃𝐹 ) − 𝑑𝑅𝑃𝐹 ) ∗ √ 𝑘𝑠𝑜𝑖𝑙 Equation 4.21 Where dRPF = depth of insulation (cm); CFDD= cumulative freezing degree day in design year (oC-day); and All other parameters are the same as before. 4.4.1.1 Example Calculate the thickness of the insulation for the given data. 1. CFDD in design year = 800 oC-day; 2. dRPF = 36 cm; 3. 𝑘𝑠𝑜𝑖𝑙 = 2.42 (W/(m.oC)); 4. The effective insulation R-value=1.0  k RPF =1/(1*100)= 0.01 (W/(m.oC)); The thickness of the insulation layer can be calculated using Equation 4.21. 129 0.01 𝑡𝑅𝑃𝐹 = (((−0.7392 ∗ 2.42 + 6.4408) ∗ (800)(.0035∗2.42+0.4846) ) − 36) ∗ √ = 5.7 𝑐𝑚 2.42 4.4.2 Gilpin’s Frost Heave Model Different theories and models for modeling frost heave are reviewed in Chapter 2 of this report. In this study, the Gilpin’s model, which is based on frozen fringe theory, was used to predict the frost heave under field conditions. As stated before the original Gilpin model is a mechanistic-empirical model based on heat and mass balance equations and laboratory data. The Gilpin model is a laboratory based model; applying it to field conditions having different boundary values led to some errors in the results. Further, the required input data to the model are not available and are expensive to obtain. Therefore, in this study, the model was simplified to include the empirical frost depth prediction model developed in this study. The resulting model was verified by comparing the predicted frost heave under pavements and shoulders to the measured values at 5 different sites in Oakland County, Michigan. 4.4.2.1 Basic Assumptions Assume that a saturated and salt-free soil column was subjected to a constant overburden pressure (POB) as shown in Figure 4.32. The top of the column was subjected to a fixed subfreezing temperature (TTOP), whereas the bottom of the column (at the ground water table elevation) was at a fixed above freezing temperature (TBOT). The soil column was further assumed to consist of three zones; frozen zone at the top followed by a frozen fringe zone and then by an unfrozen zone. The top of the unfrozen zone begins at a point where water and ice can exist in the pore spaces of the soil at below freezing temperature (Tf). In this model, frost penetration and frost heave were predicted using analytical iteration solution. In each iteration, it was assumed that 130 1. The temperature variation in each zone is linear. 2. The thermal conductivity in each zone is constant. 3. The water content and permeability in the unfrozen zone is constant. 4. Steady state water flows through the frozen fringe and unfrozen zones. POB X TTOP Frozen Zone VH H Bottom of ‘Active’ Ice Lens Tl Frozen Fringe a Tff Further Advance of Pre Ice Unfrozen Zone Z TBOT PL=0 Figure 4.32 The schematic of frost heave model (Gilpin, 1980) 131 4.4.2.2 Heat and Mass Balance Equations As stated in Chapter 2, for simulating the heat transfer in his model, Gilpin used the phasechange heat transfer equations. After imposing the boundary conditions in each zone, Equation 2-42 and 2-43 were obtained for heat transfer between the frozen and the frozen fringe zones and between the frozen fringe and the unfrozen zones, respectively. For convenience, these equations were converted to English system as follow 𝑘𝑓 (𝑇𝑇𝑂𝑃 − 𝑇𝑙 ) 𝑘𝑓𝑓 (𝑇𝑓𝑓 − 𝑇𝑙 ) 𝐿 − = 𝑉𝐻 𝐻 𝑎 𝜈𝑖 Equation 4.22 𝑘𝑓𝑓 (𝑇𝑓𝑓 − 𝑇𝑙 ) 𝑘𝑢𝑓 (𝑇𝐵𝑂𝑇 − 𝑇𝑓 ) 𝑑𝑧 − = 𝜌𝑠𝑖 𝐿 𝑎 𝑍 𝑑𝑡 Equation 4.23 − Where a = thickness of the frozen fringe (m); H = thickness of frozen zone (m); kf = thermal conductivity of the frozen zone (W/(oC.m)); kff = thermal conductivity of the partially frozen zone (W/(oC.m)); kuf = thermal conductivity of the unfrozen zone (W/(oC.m)); L = latent heat of fusion of water (J/Kg); TTOP, TBOT = temperatures at the top and bottom of the soil column (oC); VH = frost heave rate (m/s); Z = distance between bottom of soil column and position of ice penetration (m); Tl = temperature at the base of the active ice lens (oC); Tff = temperature at the base of the frozen fringe (oC); νi = specific volumes of ice (m3/kg); 𝑑𝑧 𝑑𝑡 = frost depth propagation rate (m/s); and 132 𝜌𝑠𝑖 = mass of ice per unit volume of soil (kg/m3). Using the mass balance equation and imposing boundary conditions, Gilpin proposed Equation 2-45 for calculating the water pressure at the bottom of the frozen fringe zone. Finally, Equation 2-46 was obtained for frost heave calculation. For convenience, the equations were converted to English system as follows: 𝑃𝑤𝑓 𝑉𝐻 = 𝑑𝑧 𝑍 𝜈𝑤 (𝑉𝐻 + 𝜌𝑠𝑖 Δ𝜈 𝑑𝑡 ) = −𝑔 (1 + ) 𝜈𝑤 𝜈𝑖 𝐾𝑢𝑓 𝜈𝑖2 𝑔𝜈𝑤 𝐿(−𝑇𝑙 ) [ − 𝑃𝑂𝐵 + 𝑃𝐿𝑓 ] 𝑎𝐼𝑓𝑙 𝜈𝑤 𝑇𝑚 1 [𝑇 − 𝑇 ] + (𝐾 ) 𝑓𝑓 𝑙 𝐿 1 Equation 4.24 Equation 4.25 𝑇𝑙 1 𝑑𝑇 𝑇𝑓𝑓 𝐾𝑓𝑓 𝐼𝑓𝑙 = ∫ Equation 4.26 Where 𝑔 = acceleration of gravity (m/s2); 𝑃𝑤𝑓 = water pressure at the edge of the frozen fringe (Pa); νw = specific volumes of water (m3/kg); Δ𝜈 = specific volume difference (νi − 𝜈𝑤 ); POB = overburden pressure (Pa); KL = permeability of ice lenses (m/s); 𝐾𝑢𝑓 = permeability of unfrozen zone (m/s); and All other parameters are the same as before. It should be noted that Gilpin proposed semi-empirical models for estimating the hydraulic conductivity of frozen fringe and the temperature at the bottom of the frozen fringe zone, Tf (Gilpin 1980). 133 4.4.2.3 Ice Pressure Distribution in the Frozen Fringe Zone Gilpin calculated the ice pressure distribution in the frozen fringe zone in order to model the initiation of new ice lenses. He assumed that the initiation of new ice lenses takes place where the ice pressure in the frozen fringe zone exceeds the critical pressure, which is also known as the separation pressure. This pressure is a function of the overburden pressure and water-ice curvature. Figure 4.33 illustrates the ice pressure distribution. Based on Clapeyron equation, at zero flow rate in the frozen fringe zone, the ice pressure increases along the solid line (L(-T)/(vsTa)). Further, at non-zero flow rate the ice pressure increases along the Ps line so that it becomes equal to the overburden pressure at the top of the frozen fringe zone. Ice pressure in the frozen fringe zone could be estimated using Equation 4.27. 𝑉𝐻 ∗ 𝑣𝑤 𝜈𝑖 𝑑 𝐿𝑇 = 𝐾𝑓𝑓 [𝑃𝑖 + ] 𝑣𝑖 𝑔 𝑑𝑥 𝜈𝑖 (𝑇𝑚 + 273) Equation 4.27 Where 𝐾𝑓𝑓 = the permeability in the frozen fringe (m/s); Pi = pressure in ice (Pa); Tm = bulk freezing temperature (oC); T= temperature along the frozen fringe (oC.); and All other parameters are the same as before. For modeling the separation pressure, a pair of spherical soil particles was considered as shown in Figure 4.34. In the absence of ice, the overburden pressure is acting on the interface of the particles. This pressure could be transmitted from one particle to the other. However, at a critical ice pressure, the pressure at the contact point drops to zero allowing the ice to separate the particle from each other. This critical ice pressure was assumed to be the separation pressure. Where ice pressure in frozen fringe zone exceeds the separation pressure, new ice lenses are 134 formed. Estimation of the separation pressure was a matter of debate between researchers; Gilpin suggested the following equation for separation pressure: Figure 4.33 Ice pressure along frozen fringe zone (Gilpin, 1980) 𝑃𝑠𝑒𝑝 = 𝑃𝑂𝐵 + 2 ∗ 𝜎𝑖𝑤 𝐷10 Where 𝑃𝑠𝑒𝑝 = separation pressure (Pa); D10= = particle size at 10 percent passing (m); 𝜎𝑖𝑤 = ice-water surface energy (N/m); and 𝑃𝑂𝐵 = overburden pressure (Pa). 135 Equation 4.28 Other equations were developed and are available in the literature. However, in this study, the Gilpin equation was revised to simplify the required inputs and used to predict the frost heave potential. R Ice Pore Water Particle Figure 4.34 Particle separation pressure 4.4.3 Revised Frost Heave Model In the Gilpin model at the beginning of the solution, initial non-zero values were chosen for a and H parameters of Equation 4.22 and Equation 4.23 in order to avoid the infinite temperature gradient (See Figure 4.32). At each iteration (Δt), the systems of four equations (Equation 4.22 to Equation 4.25) were solved to calculate the four unknowns, i.e. VH, Pwf, Tl and dz/dt. It should be noted that since Equation 4.25 is nonlinear the accuracy of the results are highly related to the nonlinear solution. After calculating VH, ice pressure in the frozen fringe was calculated using Equation 4.27. If the ice pressure did not exceed the critical pressure then H was increased by VH*Δt; a was increased by dz/dt* Δt and the equations were solved for the next iteration. Otherwise a new ice lens was assumed to initiate where ice pressure exceeded the 136 critical pressure, H was increased and a was decreased accordingly. Then VH, Pwf, Tl and dz/dt were calculated again for the same time step (Gilpin, 1980). In Gilpin model, the hydraulic conductivity of the frozen fringe zone (Kf) was estimated based on the laboratory data. Since laboratory conditions were not necessarily correlated well with the field conditions and field data were not available for calibration, an overall permeability (Kf) was assumed for the frozen fringe zone in order to avoid nonlinear solution. Further, due to a large frozen zone thickness in the field, the frozen zone was assumed to be impermeable. Tff (temperature at the bottom of frozen fringe) was calculated using the following empirical equation (Gilpin 1980). 𝑇𝑓𝑓 = − 8𝜎𝑖𝑤 𝜈𝑤 𝑇𝑚 𝐷10 ∗ 𝐿 Equation 4.29 Where all parameters are the same as before. In the revised model, instead of using Equation 4.23 for calculating the frost depth propagation, the empirical frost depth model that was developed based on the measured frost depths data in Michigan was used (Equation 4.12). The analyses were conducted using analytical iterative solution. In each iteration, the frost depth propagation rate ( 𝑑𝑧 𝑑𝑡 ) was calculated using Equation 4.30 and Equation 4.31 and the Tl was estimated using Equation 4.32. 𝑑𝑧 = (−0.7392 𝑘𝑢𝑓 + 6.4408) 𝛽 𝑑𝑡 Equation 4.30 𝛽 = (𝐶𝐹𝐷𝐷𝑛 (0.035𝑘𝑢𝑓+0.4846) − 𝐶𝐹𝐷𝐷𝑛−1 (0.0035𝑘𝑢𝑓+0.4846) ) Equation 4.31 𝑇𝑙 = 𝑇𝑓𝑓 − 𝑎 𝑑𝑧 𝑘𝑢𝑓 (𝑇𝐵𝑂𝑇 − 𝑇𝑓 ) ∗ (𝜌𝑠𝑖 𝐿 + ) 𝑘𝑓𝑓 𝑑𝑡 𝑍 Where CFDD= Cumulative freezing degree day (oC- day); and All parameters are the same as before 137 Equation 4.32 Furthermore, Equation 4.25 was revised into Equation 4.33 and by solving Equation 4.24 and Equation 4.33; the values of VH and Pwf were calculated. 𝜈𝑖2 𝐾𝑓𝑓 𝐿(−𝑇𝑙 ) 𝑉𝐻 = [ − 𝑃𝑂𝐵 + 𝑃𝑤𝑓 ] 𝑔𝜈𝑤 𝑎 𝜈𝑤 (𝑇𝑚 + 273) Equation 4.33 Where 𝐾𝑓𝑓 = over all permeability of frozen fringe (m/s); and All other parameters are the same as before. Finally, the ice pressure variation in the frozen fringe zone was calculated using Equation 4.27. New ice lens formation was assumed where the pressure value is higher than the separation pressure. Therefore, the thicknesses of the frozen fringe and the frozen zone were changed accordingly and consequently, the calculations of VH, Tl, Pwf were repeated. The total frost heave was then estimated using the following equation ∆ℎ𝑡𝑜𝑡𝑎𝑙 = ∆ℎ𝑢 + ∆ℎ𝑖 = 𝑉𝐻 ∗ ∆𝑡 + 0.09 ∗ 𝑛 ∗ 𝐻 Equation 4.34 Where ∆ℎ𝑡𝑜𝑡𝑎𝑙 = total frost heave (m); ∆hu = frost heave due to water uptake (m); ∆hi = heave due to freezing of in-situ pore water (m); ∆t = time interval (s); n = soil porosity, and 0.09 is the ratio of volumetric expansion of water in phase change (Nixon 1982); and all other parameters are the same as before. 4.4.4 Discussion of the Results of the Revised Frost Heave Model At the beginning of the frost, the heave rate is high therefore the ice pressure in the frozen fringe zone could surpass the separation pressure and the boundary of the frozen zone keeps moving downward in the soil column. As frost progresses, the heave rate decreases and consequently the ice pressure decreases and it does not exceed the separation pressure anymore. 138 This causes growth in the frozen fringe zone. Since the hydraulic conductivity of the frozen fringe zone is less than that of the unfrozen zone the larger frozen fringe zone leads to lower heave. The extent of the frozen fringe zone depends on the types of the soil and overburden pressure. Soil Type - The hydraulic conductivity of fine sand is higher than the clayey silt so the flow rate is greater in fine sand. However, as mentioned before, capillary pressure is mainly responsible for the frost heave phenomenon. Due to the aggregate size, suction is smaller in fine sand than in clayey silt, therefore larger frost heave rates are expected in clayey silt. Further, in fine sands relative to clayey silts a larger frozen fringe zone is observed. Figure 4.35 and Figure 4.36 depicts the calculated frost heave and frozen fringe thickness in three different types of soil versus time, respectively. Clayey silt Sandy clayey silt Fine sand 35 30 Frost heave (mm) 1. 25 20 15 10 5 0 10 30 50 Time (Day) 70 90 Figure 4.35 Calculated frost heave for three soil types when GWL= 10 m, TTOP= -3 oC for 100 days, POB = 150 kPa 139 Clayey silt 120 Sandy clayey silt Fine sand Frozen fringe(mm) 100 80 60 40 20 0 10 30 50 Time (Day) 70 90 Figure 4.36 Calculated frozen fringe for three soil types when Z= 10 m, TTOP=- 3 oC for 100 days, POB = 150 kPa 2. Overburden Pressure - Various overburden pressures were used to assess the impact of overburden pressure on the frost heave and frozen fringe zone. Konrad and Morgenstern (1982) found out that the overall permeability of the frozen fringe zone decreases approximately by 25% as the overburden pressure increases up to 400kPa. Accordingly, in the revised model, the frozen fringe zone permeability was reduced as the overburden pressure was increased. Figure 4.37 to Figure 4.40 show the model results for total heave in different overburden pressures for different Z values when the TTOP was fixed at -3 oC for 100 days in different soil types. As can be seen in the figures, when the ground water table is deep, for the same freezing time period, less amount of water can reach the frozen zone and therefore the frost heave decreases. 140 6 z=6 m z= 7.5 m z=9 m Total heave (cm) 5 4 3 2 1 0 0 50 100 150 200 Overburden pressure (kPa) 250 300 Figure 4.37 Calculated total heave versus overburden pressure for clayey silt in different ground water table depths when TTOP=-3 oC in 100 days 5 z=6 m z= 7.5 m z=9 m Total heave (cm) 4 3 2 1 0 0 50 100 150 Overburden pressure (kPa) 200 Figure 4.38 Calculated total heave versus overburden pressure for sandy clayey silt in different ground water table depths when TTOP=-3 oC in 100 days 141 4.5 z=6 m z= 7.5 m z=9 m 4 Total heave (cm) 3.5 3 2.5 2 1.5 1 0.5 0 0 50 100 150 Overburden Pressure (kPa) 200 Figure 4.39 Calculated total heave versus overburden pressure for fine sand and silt with pebbles in different ground water table depths when TTOP=-3 oC in 100 days 4.5 z=6 m z= 7.5 m z=9 m 4 Total heave (cm) 3.5 3 2.5 2 1.5 1 0.5 0 0 50 100 Overburden Pressure (kPa) 150 Figure 4.40 Calculated total heave versus overburden pressure for clayey, silty, gravely, sand in different ground water table depths when TTOP=-3 oC in 100 days 142 At higher overburden pressures, the separation pressure is larger. Therefore, it is expected that the frozen fringe zone thickness increases as the overburden pressure increases. This leads to lower frost heave values in higher overburden pressures. In order to assess the impact of surface temperature (i.e. TTOP), the model was assessed with a fixed TTOP and with a changing TTOP but a fixed rate of cooling. In both scenarios, the freezing index was the same. The results are depicted in Figure 4.41. It can be seen from the figure that the results are almost the same in both assessments, therefore using a fixed TTOP based on the cumulative freezing index and length of frost period is a good assumption. 7 Variable Ttop Fixed Ttop Total heav (cm) 6 5 4 3 2 1 0 0 100 200 300 Overburden pressure (Kpa) 400 Figure 4.41 Total heave versus overburden pressure for clayey silt when TTOP is fixed at -3 oC and when TTOP is decreasing with a rate of -.057 per day in 100 days Furthermore, the revised model was used to predict the frost heave under shoulder and pavement. Data from 5 sites located in Oakland County, Michigan were provided by the 143 Michigan Department of Transportation (See Chapter 3-4). In order to validate the revised model, this data was used. According to each soil description and size distribution, hydraulic conductivity and D10 values were chosen. Thermal conductivity values were chosen according to the soil type from the measured thermal conductivity values (See Table 3-4). Furthermore, by using the measured frost depth at maximum heave and Equation 4.12, the CFDD values and the TTOP for each site were calculated. The inputs for each soil type are shown in Table 4.10. As stated before, the frost heave occurs when water migrates from the water table to the frozen layer. Therefore, the controlling layer is the natural soil. It can be assumed that frost heave could be calculated by a single layer model consisting of the natural soil under the pavement layers. The weight of the asphalt, base, and subbase layers can be considered as the overburden pressure. The frost heave was estimated under the shoulder and pavement in Oakland County sites. The results are shown in Figure 4.42. The only difference between the two models is the overburden pressure. The overburden pressure was modeled by a 51-centimeter thick soil layer having density of 19.6 (kN/m3) for the shoulder and a 77.5-centimeter thick soil layer with the same density for the pavement. It should be noted that the ground water table was set at the average measured ground water table level in Oakland County, which was 9 meters. Figure 4.42 indicates that in both cases the calculated frost heave values are within 0.25 centimeters of the measured ones. It can be concluded that different simplifications and modifications, which were applied to the Gilpin model, did not affect the accuracy of the model significantly, indeed, it produced better results. 144 Table 4.10 Different input values for each site, I75, Oakland County, Michigan. Station Name Duration (days) Ttop (oC) Sta/724+00 65 -1.94 Hydraulic Conductivity (m/s) 1.0*10^-6 Sta/719+00 40 -3.06 Sta/652+00 60 Sta/528+88 Sta/474+00 9 D10 (mm) 0.01 1.0*10^-6 9 0.01 -3.89 5.0*10^-7 9 0.02 70 -1.94 1.0*10^-7 9 0.002 55 -1.94 5.0*10^-8 9 0.001 Line of equality between calculated and measured data Frost heave under shoulder Frost heave under pavement 3 Calculated heave (cm) GWT (m) 2.5 2 1.5 1 0.5 0.5 1 1.5 2 Measured heave (cm) 2.5 3 Figure 4.42 Measured versus calculated frost heave under the shoulder and pavement in 5 sites, Oakland County, Michigan 145 It is noteworthy that for station 528+88 (See Table 4.10); the measured frost heave under the pavement is approximately 0.9 centimeters less than that under the shoulder (See Table 3-5). This difference is higher than those at the other sites (about 0.3-centimeter). At station 528+88, an undercut of approximately 30.5 centimeters was made for frost protection of the roadbed soil. Therefore, the frost penetration in the clayey silt roadbed soil decreased by 30.5 centimeters and hence, as it was expected, the frost heave decreased. In the analyses, the undercut was modeled as a part of the overburden pressure against the surface of the clayey silt roadbed soil. The results are also shown in Figure 4.42. 4.4.5 Heave Pressure Heave pressure can cause real stability issues in different structures such as retaining walls, utility poles, and shallow foundations. In the revised frost heave model presented in the previous section, frost heave can be calculated as a function of overburden pressure. If the overburden pressure is equal to or greater than the heave pressure, no heave will occur. That is equilibrium scenario is reached. Otherwise, frost heave will take place due to the net pressure against the structure in question. The heave pressure can be calculated as follows: 𝑃𝐹𝐻 = 𝑃𝐸 − 𝑃𝑂𝐵 Equation 4.35 Where PFH = pressure due to heave (kPa); PE = the equilibrium overburden pressure (see Figure 4.37), (kPa); POB= the actual overburden pressure (kPa); In order to develop a model for estimating the frost heave pressure, the following three steps were used: 146 1. For four soil types, the revised heave model (Equation 4.30 to Equation 4.34) was used to calculate the amount of heave as a function of the overburden pressure and the depth to the ground water table (see Figure 4.37 through Figure 4.39). Table 4.11 shows the different input values for each soil type. For each scenario, the corresponding equilibrium pressure was also calculated. Table 4.11 Different input values for each soil type Soil Type Hydraulic D10 Conductivity (mm) (Kuf, (m/s) Thermal Conductivity (W/(m.oC)) Dry unit Weight (γd, kN/m3) Water Content (w%) POB (kPa) Clayey, silty, 1.0*10^-6 0.02 2.57 10 19.6 gravely, sand Fine sand variable 5.0*10^-7 0.01 2.42 15 18.9 and silt Sandy 1.0*10^-7 0.002 1.75 20 18.1 clayey silt clayey silt 5.0*10^-8 0.001 1.52 25 15.7 1- Soil type: can be obtained from the boring log on site (known) 2- Duration: the time period that CFDD (cumulative freezing degree day) is calculated over (known or assumed) 3- Ttop: temperature at the ground surface= CFDD/Duration (known or assumed) 4- Hydraulic conductivity: can be measured on site or assumed based on the soil type 5- GWTD: ground water table depth (known) 6- D10: the effective size of the soil; can be obtained from the soil distribution curve or assumed based on the soil type (known or assumed) 7- Thermal conductivity: measured at MSU soil laboratory (known) 8- Tbottom: temperature at the ground water table level; assumed based on GWTD (assumed) 9- Dry unit weight of soil: can be measured on site or obtained from the CRREL graphs based on the thermal conductivity values (known or assumed) 10- Water content: can be measured on site or obtained from the CRREL graphs based on the thermal conductivity values (known or assumed) 11- Void ratio: can be measured on laboratory or assumed based on the soil type and its density= 0.5 (known or assumed) 12- POB= overburden pressure 147 2. The data in Figure 4.37 through Figure 4.39 were used to estimate (for each amount of heave and ground water depth), the corresponding overburden pressure. 3. The estimated overburden pressure and the equilibrium pressure were used as inputs to Equation 4.35 to estimate the heave pressure. Figure 4.43 shows the results for clayey silt. As can be seen the heave pressure is almost the same in different ground water table depth. Therefore, heave pressure can be estimated regardless of the ground water table depth as a function of frost heave. Figure 4.44 shows the heave pressure versus frost heave in four soil types. 6 z=6 m z=7.5 m z=9 m Total heave (cm) 5 4 3 2 1 0 0 50 100 150 200 250 300 Frost pressure (kPa) Figure 4.43 heave pressure versus calculated total heave for clayey silt in different ground water table depths when TTOP= -3 oC in 100 days 148 The results showed that in the same winter duration the heave pressure has a unique polynomial relationship with frost heave as follows 2 𝑃𝐹𝐻 = 𝑎 ∗ ∆ℎ𝑡𝑜𝑡𝑎𝑙 + 𝑏 ∗ ∆ℎ𝑡𝑜𝑡𝑎𝑙 Equation 4.36 Where a,b,c = constant values which are different in each soil type; and all other parameters are the same. It should be noted that in Equation 4.36 the soils were considered to be saturated. Also, the effect of void ratio was not considered in the model (since the void ratio is a function of soil density) Table 4.12 shows the statistical parameters of Equation 4.36 for each soil type. Clayey Silt Sandy Clayey Silt Clayey, Silty, Gravely Sand Fine Sand with Silt 250 Heave pressure (kPa) 200 150 100 50 0 0 1 2 3 Total heave (cm) 4 5 6 Figure 4.44 Heave pressure versus calculated total heave in four soil types when TTOP= -3 oC over 100 days 149 Table 4.12 Statistical coefficients in Equation 4-30 for each soil type. Soil Type Equation 4.36 clayey silt PFH = -2.86∆htotal2 + 60.2∆htotal Sandy clayey silt PFH = -2.64∆htotal 2 + 49.1∆htotal Fine sand and silt PFH = -2.63∆htotal 2 + 45.4∆htotal clayey, silty, gravely, sand PFH = -2.17∆htotal 2 + 37.1∆htotal It should be noted that field data for evaluating the accuracy of the results were not available. Therefore, the model should be validated as data become available. Thaw Depth At the end of the freezing season, the soils start to thaw. The prediction of frost and thaw depths are crucial for estimating the amount of heave due to frost action and to estimate the proper time to post and remove seasonal load restriction signs. The calculation of thaw depth is presented and discussed below. 4.5.1 Calculation of Cumulative Thawing Degree day (CTDD) Calculating the CTDD accurately is the first step in developing an accurate thaw depth model and consequently an effective and reliable SLR policy. In this study, the average daily air and the pavement surface temperatures in 2011(which were available at the early stage of the study) from 12 Road Weather Information System (RWIS) stations in Lower Peninsula (LP) of the State of Michigan were used to develop a new CTDD approach. First, surface and air temperature data and the calculated solar radiation were used to develop a model for estimating the average daily pavement surface temperature. Second, the estimated surface temperatures were used to calculate the CTDD. Later on, when more data became available, the 2012, 2014, and 2015 data from the same RWIS stations were used to check the accuracy of the predicted surface temperatures for different years with different winter severity. Third, the accuracy of the model was further evaluated using the air and the average daily pavement surface temperatures 150 in 2012, 2014, and 2015 from 9 RWIS stations in the Upper Peninsula (UP) of the State of Michigan. The results of the model were also compared with WSDOT and MnDOT predictions (See Chapter 2, Section 2.3.2). 4.5.1.1 Pavement Surface Modeling For modeling the daily pavement surface temperature, the air-pavement system can be considered as a thermodynamic system and the energy balance (equilibrium) at the pavement surface can be calculated using Equation 2.2. As mentioned in Chapter 2, Solaimanian and Bolzan proposed a mechanistic empirical nonlinear model for solving Equation 2.2 (Solaimanian and Bolzan, 1993). The model requires various inputs such as air temperature, solar radiation, wind speed, percent of sunshine, pavement emissivity, pavement absorptivity, and so forth. Since most of these input values are not available and/or expensive to collect, most highway agencies only consider air temperature as an input in their calculations of the pavement surface models and CTDD. However, Solaimanian and Bolzan results showed that while the air temperature significantly affects the pavement surface temperature predictions, the difference between air temperature and pavement temperature can be as low as 5.5 to 8 oC or as high as 22 to 28 oC depending on the solar radiation and percent sunshine. The solar radiation can change based on the latitude and time of the year. In fact, different locations of earth receive different amount of solar radiation due to solar inclination, i.e. the tilt of the earth north-south axis with respect to the orbital plane. On the other hand, during the year the distance between the earth and the sun changes, which causes a daily variation in the amount of solar radiation (Diefenderfer et al, 2006). Therefore, in this study, in addition to air temperature, calculated solar radiation was also considered as an input of the developed model. At any geographical location, the daily amount of solar radiation was 151 calculated using the latitude and the day of the year as follows (Diefenderfer et al, 2006; Iqbal, 2012) 𝜔𝑠 ∗ 𝜋 − tan 𝜔𝑠 ) 180𝑜 𝐷𝐶𝑆𝑅 = (76.39) ∗ 𝐼𝑠𝑐 ∗ 𝐸0 ∗ sin(𝜑) ∗ sin(𝛿) ∗ ( Equation 4.37 𝑑𝑛 ] 365 Equation 4.38 𝐸0 = 1 + 0.033 ∗ cos [2𝜋 𝛿 = 23.45𝑜 ∗ sin [ 360𝑜 ∗ (𝑑𝑛 + 284)] 365 𝜔𝑠 = cos−1 (− tan(𝜑) ∗ tan(𝛿)) Equation 4.39 Equation 4.40 Where DCSR = daily calculated solar radiation ((kJ/(m2day)); Isc = solar constant= 4,871 (kJ/(m2.hr)) E0 = daily eccentricity factor; φ = latitude (deg); δ = solar declination (deg); ωs= sunrise hour angle which is the angle between the sun’s highest point each day, which is zero, and the location of the sun at sunrise in that day(deg); and dn = day of the year starting at 1 for January first; for example, February 4 is the 35th day of the year, hence, the corresponding dn is 35. In order to develop a prediction model for pavement surface temperature, the latitude of each RWIS station and the day of the year (dn) were used as inputs to Equation 4.37 to calculate the solar radiation for each day where the air and pavement surface temperatures were measured. Two statistical pavement surface temperature models were then developed; linear and nonlinear as stated in Equation 4.41 and Equation 4.42. For both models, the measured pavement surface temperature was treated as a dependent variable and the output of Equation 4.37 (the solar 152 radiation), and the average measured daily air temperature as the independent variables. The metrics of the two models are depicted in Table 4.13 𝑇𝑠𝑢𝑟𝑓 = 0.687𝑇𝑎𝑖𝑟 + 4.377 ∗ 10−4 ∗ 𝐷𝐶𝑆𝑅 − 6.4 Equation 4.41 𝑇𝑠𝑢𝑟𝑓 = 0.674𝑇𝑎𝑖𝑟 + 1.61 ∗ 10−4 ∗ 𝐷𝐶𝑆𝑅1.09 − 5.8 Equation 4.42 Where Tsurf = average daily pavement surface temperature (oC); Tair = average daily air temperature (oC); and all other parameters are the same as before. Table 4.13 Performance metrics for the model Model R2 Standard Error of the Estimate (SEE) Mean Absolute Error(MAE) Equation 4.41 0.89 3.6 2.9 Equation 4.42 0.91 2.8 2.2 Figure 4.45 depicts the measured and the calculated pavement surface temperatures using Equation 4.41 and Equation 4.42 for 12 RWIS stations in 2011. It should be noted that since in most years the pavement is completely thawed in the first 120 days of the year, in this figure data from 12 RWIS stations of the first 120 days of the year in 10-day interval are shown. The solid straight line in Figure 4.45 indicates the locus of equality between the measured and the calculated pavement surface temperatures. It can be seen from the figure that Equation 4.42 predicts the pavement surface temperature slightly better than Equation 4.41 especially for pavement surface temperature higher than 15.5oC. Therefore, it is recommended to use Equation 4.41 for more accurate prediction of the pavement surface temperature. 153 Equation 4.41 Equation 4.42 Line of Equality Calculated surface temperature (oC) 25 20 15 10 5 0 -15 -10 -5 0 5 10 15 20 25 -5 -10 -15 Measured surface temperature (oC) Figure 4.45 Calculated versus measured surface temperature in 2011 using Equation 4.41 and Equation 4.42 for 12 RWIS stations in the first 120 days of the year. The data are shown in 10day interval After developing the statistical model for predicting the pavement surface temperature, Equation 4.42 was used to calculate, for each day, the thawing degree-day using Equation 4.43. Finally, the cumulative thawing index was calculated using Equation 4.44. 𝑇ℎ𝑎𝑤𝑖𝑛𝑔 𝐷𝑒𝑔𝑟𝑒𝑒 𝑑𝑎𝑦 = (𝑇𝑠𝑢𝑟𝑓 − 0𝑜 𝐶) 𝑛 Equation 4.43 Equation 4.44 𝐶𝑇𝐷𝐷𝑛 = ∑ 𝑇ℎ𝑎𝑤𝑖𝑛𝑔 𝐷𝑒𝑔𝑟𝑒𝑒 𝐷𝑎𝑦 ≥ 0 𝑖=1 Where 𝐶𝑇𝐷𝐷𝑛 = Cumulative thawing Index inform January first to the nth day of the year (°Cday); and All other parameters are the same as before. 154 4.5.1.2 Model Verification In order to evaluate the accuracy of the statistical model, Equation 4.42 was used to calculate the pavement surface temperatures for three years other than the year 2011 because the 2011 data were used to develop the statistical model. The three years (2012, 2014, and 2015) have substantially different winter characteristics and severity. For example, the 2015 winter was an average winter, 2012 winter was a relatively warm one, and 2014 winter was a particularly cold winter (Based on the National Oceanic and Atmospheric Administration (NOAA) weather data, the Normal freezing index for LP and UP is 360 and 555 degree-days, respectively. Therefore, the year with approximately the same freezing index counted as average, the year with the lower freezing index counted as warm and the year with the higher freezing index counted as cold). Nevertheless, Equation 4.42 was used to calculate the average daily pavement surface temperature for the three years and for the 12 RWIS stations in the LP of the State of Michigan. The latitudes of these stations are approximately in the range of 43.9o to 45.8o. The results are depicted in Figure 4.46. In the figure, the number of points (n) and the coefficient of determination (R2) are shown. The data in the figure indicates that Equation 4.42 accurately predicts the measured pavement surface temperature for all three years and for the 12 RWIS in the Lower Peninsula of the State of Michigan. Indeed, the coefficient of determination (R2) of Equation 4.42 is slightly higher than the value of R2 for the 2011 data, which were used in developing Equation 4.42. In addition, Equation 4.42 was also used to predict the measured pavement surface temperatures at 9 additional RWIS stations located in the Upper Peninsula (UP) of the State of Michigan. The latitudes of these stations are approximately in the range of 45.9o to 46.8o. The results are shown in Figure 4.47. It can be seen from the figure that Equation 4.42 accurately predicted the daily pavement surface temperatures for most days. 155 Line of Equality Calculated surface temperature(oC) 30 -20 R2=.94 n=3294 20 10 0 -10 0 10 20 30 -10 -20 Measured surface temperature(oC) Figure 4.46 Calculated versus measured surface temperature in 2012, 2014 and 2015 for 12 RWIS stations in the LP of the State of Michigan Calculated Surface Temperature(oC) Line of Equality -20 30 R2=0.91 n=1691 20 10 0 -10 0 10 20 30 -10 -20 Measured Surface Temperature(oC) Figure 4.47 Calculated versus measured surface temperature in 2012, 2014 and 2015 for 9 RWIS stations in the UP of the State of Michigan 156 Furthermore, Equation 4.44 was used to calculate the CTDD for the measured and the calculated pavement surface temperature. Figure 4.48 depicts the results for the first 120 days of the year in 2012, 2014 and 2015 for 12 RWIS stations in the LP of the state of Michigan. In this figure, the number of data points (n) and the mean absolute percentage error (MAPE) are shown. In addition, Equation 4.44 was also used to calculate the CTDD for the measured and the calculated pavement surface temperature in 2012, 2014 and 2015 for 9 RWIS stations in the UP of the State of Michigan. The results are shown in Figure 4.49. It can be seen that Equation 4.44 predicts the CTDD less accurately for the UP stations. In fact, the MAPE of Equation 4.44 is 18% for the LP stations and 23% for the UP stations. Additionally, the WsDOT and MnDOT methods were used to calculate the average daily pavement surface temperature in 2012, 2014, and 2015 for 12 RWIS stations in the LP and 9 RWIS stations in the UP of the state of Michigan. Table 4.14 shows the different performance metrics for each model. In this table, the standard error of the estimate (SEE), the mean absolute error (MAE) and the mean absolute percentage error (MAPE) are shown for each model. As can be seen in Table 4.14, for the LP stations, the SEE values of MnDOT and WSDOT models are about 28% and 60% larger than SEE of Equation 4.42, respectively. This difference is even larger for the UP stations. This indicates that Equation 4.42 predicts the average pavement surface temperature better than the other two methods. In general, all of the performance metrics indicate that Equation 4.42 predicts the average pavement surface temperature with higher accuracy. Further, the accuracy of the calculated CTDD was also compared to the accuracy of the WSDOT and MnDOT methods. The CTDD values for the years 2012, 2014, and 2015 and for 157 each of the 12 RWIS stations in the LP and 9 RWIS stations in the UP were calculated using the three methods. Table 4.15 shows various performance metrics for each model. Line of Equality CTDD of calculated surface temperature (oC-day) 700 MAPE= 18% n-= 1677 600 500 400 300 200 100 0 0 100 200 300 400 500 600 o CTDD of measured temperature( C-day) 700 Figure 4.48 CTDD of the calculated surface temperature versus CTDD of the measured surface temperature in 2012, 2014, and 2015 for 12 RWIS stations in the LP of the State of Michigan CTDD of calculated surface temperature (oC-day) Line of Equality 600 MAPE= 23% n= 718 500 400 300 200 100 0 0 100 200 300 400 500 600 o CTDD of measured surface temperature ( C-day) Figure 4.49 CTDD of the calculated surface temperature versus CTDD of the measured surface temperature in 2012, 2014, and 2015 for 9 RWIS stations in the UP of the State of Michigan 158 Table 4.14 Performance metrics of Equation 4.42, MNDOT and WSDOT models for the average daily pavement surface temperature for 12 RWIS stations in the LP and 9 RWIS stations in the UP in 2012, 2014, and 2015 Location Model Equation 4.42 Lower MnDOT Model Peninsula WSDOT Model Equation 4.42 Upper MnDOT Model Peninsula WSDOT Model Standard Error of the Estimate (SEE) Mean Absolute Error (MAE) 2.3 2.7 3.6 2.9 3.8 4.7 1.8 2.2 2.5 2.3 2.8 3.7 Mean Absolute Percentage Error (MAPE) 0.13 0.17 0.19 0.10 0.29 0.35 As can be seen, for all of the UP and LP stations, the standard error of the estimate (SEE) values of MnDOT and WSDOT models are, respectively, about 15% and 135% higher than the SEE of Equation 4.44. Such differences are significant and adversely affect the timely posting and removing the SLR signs. Such timing could be off by few days to about a week when the pavement is in critical condition states. Accurate prediction of posting and removing the SLR signs saves the road owners and the trucking industries unnecessary expenses. Table 4.15 Performance metrics of Equation 4.44, MnDOT and WSDOT models for the CTDD calculation for 12 RWIS stations in the LP and 9 RWIS stations in the UP in 2012, 2014, and 2015 Location Lower Peninsula Upper Peninsula Model Standard Error of the Estimate (SEE) Mean Absolute Error (MAE) Equation 4.44 MnDOT Model WSDOT Model Equation 4.44 MnDOT Model WSDOT Model 19 22 44 22 26 51 11 13 24 11 13 25 159 Mean Absolute Percentage Error (MAPE) 0.17 0.26 0.34 0.23 0.32 0.48 4.5.2 Nixon and McRoberts Equation After calculating the CTDD, Nixon and McRoberts equation (Equation 2-87) was used to estimate the depth of thaw at the various RWIS stations in the state of Michigan. Figure 4.50 and Table 4.16 depict the results. It can be seen from the figure that the results are not satisfactory. In fact, Equation 2-87 under predicts thaw depth by as much as 76 centimeters in some stations. The error could be related to the simplifying assumptions made in the equation, the lack of exact input data, or error in calculating the thaw index. Since Nixon and McRoberts Equation did not yield accurate thaw depth results, new empirical models were developed in this study using the RWIS data in the State of Michigan. Line of equality between the measured and calculated data Maximum calculated thaw depth (cm) 160 140 120 100 80 60 40 20 0 0 20 40 60 80 100 120 140 Maximum measured thaw depth (cm) 160 Figure 4.50 Maximum thaw depths predicted by Nixon and McRoberts equation versus the measured maximum thaw depths in Michigan. 160 Table 4.16 Maximum thaw depth predicted by Nixon and McRoberts equation for RWIS stations Location Lower Peninsula Name of the Station Benzonia Cadillac Grayling Houghton Lake Reed City Waters Au Train Brevort Upper Peninsula Harvey Golden Lake Seney Twin Lakes Type of Soil Year Maximum Measured Thaw Depth (cm) Loose Sand Dense Sand Dense Sand 2010-2011 2010-2011 2010-2011 53 53 117 26 25 50 Dense Sand 2010-2011 117 62 2010-2011 71 24 2010-2011 53 34 2010-2011 117 48 2010-2011 117 53 2010-2011 117 56 2010-2011 157 61 2010-2011 2010-2011 102 102 58 48 Compacted Sand Loose Sand Compacted Sand Loose Sand Sand with Gravel and Silt Loose Sand Loose Sand Sand with Gravel and Silt Dense Sand Dense Sand with Gravel Loose Sand Silty Clayey Sand Nixon Maximum Calculated (cm) 4.5.3 Thaw Depth Empirical Models Since Nixon and McRoberts Equation did not yield accurate thaw depth results, new empirical models were developed in this study using the RWIS data in the State of Michigan and Minnesota. These new models are presented below. Among the 25 RWIS stations located in the Upper and Lower Peninsulas of Michigan, only 12 stations were used for developing the empirical models. Thirteen stations were not considered due to incomplete data (some of the data are missing). For air temperature data the 161 nearest NOAA station database was used due to higher consistency and accuracy with respect to the RWIS database. Sites not situated close to NOAA stations are not considered in this study. Due to data availability, data from the spring of 2011, 2014, and 2015 were used. Unfortunately, all of the 12 stations contain sandy soils. First, the air temperature data were used to calculate the CTDD for each RWIS station location. Second, the measured thaw depth data in all RWIS stations and the calculated CTDD values were used to develop a statistical prediction equation for sandy soil. This resulted in Equation 4.45 𝑇 = −0.0028 ∗ (𝐶𝑇𝐷𝐷)2 + 1.6999 ∗ 𝐶𝑇𝐷𝐷 Equation 4.45 Where T = thaw depth (cm); and CTDD = cumulative thawing degree day (oC – day). The results of Equation 4.45 are depicted in Figure 4.51. The calculated coefficient of determination (R2) is 0.85. As expected, the accuracy of the thaw depth model is lower than the frost depth model. The main reason is that the variability of the measured thaw depth data is higher than that of the frost depth data. Such variability is a function of the measured solar radiation, percent sunshine, absorptivity, emissivity, and so forth. Unfortunately, such data are not available at this time to improve Equation 4.45. Nevertheless, the equation does predict the thaw depth data relatively accurately. Figure 4.52 shows the calculated thaw depths using Equation 4.45 versus the measured ones. The solid straight line in the figure is the locus of equality between the measured and calculated data. As can be seen, the majority of the calculated thaw depth data are within a few centimeters from the measured values. 162 220 T = -0.0028CTDD2 + 1.6999CTDD R² = 0.8519 n= 143 200 Thaw depth (cm) 180 160 140 120 100 80 60 40 20 0 0 20 40 60 80 100 120 o Cumulative thawing degree day ( C-day) 140 160 Figure 4.51 Measured thaw depths versus cumulative thawing degree-day for sandy soils in the State of Michigan Calculated thaw depth (cm) Line of equality between the measured and calculated data Sandy soil-Michigan 200 180 160 140 120 100 80 60 40 20 0 n= 143 0 20 40 60 80 100 120 140 Measured thaw depth (cm) 160 180 200 Figure 4.52 Measured thaw depth in sandy soils in the state of Michigan versus the calculated values using Equation 4.45 163 As stated in Chapter 3, thaw depth data measured from 2003 to 2012 in 8 stations in the State of Minnesota were requested and received from MnDOT. Equation 4.45 was then used to calculate the thaw depth data in sandy stations and for the ten-year period. The measured thaw depth data and the calculated ones are depicted in Figure 4.53 . The straight line in the figure indicates the line of equality between the measured and the calculated values. Calculated thaw depth (cm) Line of equality between the measured and calculated data Sandy soil- Minnesota 220 200 180 160 140 120 100 80 60 40 20 0 n= 124 0 20 40 60 80 100 120 140 160 180 200 220 Measured thaw depth (cm) Figure 4.53 Measured thaw depths in sandy soil in the state of Minnesota versus the thaw depth values calculated using Equation 4.45. It should be noted that for sandy soil (Figure 4.53) the number of measured data points is 130. The calculated coefficient of determination (R2) for Minnesota data is 0.82 which is slightly lower than R2 for Michigan data. The relatively high values of R2 indicate that Equation 4.45 predict the thaw depth data in sandy soils in the states of Michigan and Minnesota relatively accurately. In order to further evaluate the performance of Equation 4.45, various performance metrics for both states are shown in Table 4.17. 164 Table 4.17 Performance metrics of Equation 4.45 for thaw depth estimations in the State of Michigan and Minnesota Location Standard Error of the Estimate (SEE) Mean Absolute Error (MAE) Michigan 20.3 16.0 Mean Absolute Percentage Error (MAPE) 0.30 Minnesota 0.29 20.6 16.3 In this table, the standard error of the estimate (SEE), the mean absolute error (MAE) and the mean absolute percentage error (MAPE) are shown for each model. The SEE is proportional to the width of the confidence interval so smaller SEE is an indication of a better fit. As can be seen in, Table 4.17 the model performance is almost identical in both states. To further evaluate the validity of Equation 4.45, a statistical model was developed using the measured thaw depth and the calculated CTDD for sandy soils in Minnesota. The results are Thaw depth (cm) shown in Figure 4.54. Statistical polynomial function Poposed model (Eq.4-45) 220 200 180 160 140 120 100 80 60 40 20 0 n=124 0 20 40 60 80 100 120 140 160 180 200 220 Cumulative thawing degree days (oC-day) Figure 4.54 thaw depths versus cumulative thawing degree-day for sandy soil showing the best fit and proposed model (Equation 4.45) in the State of Minnesota. 165 The dashed and solid curves in Figure 4.54 represent the statistical model and Equation 4.45, respectively. Examination of Figure 4.54 indicates that up to 150 centimeters of thaw depth, Equation 4.45 and the statistical model produce almost the same results for sandy soils in Minnesota. After 150 centimeters, the differences between the results of the two models are less than 20 centimeters. In any events, the number of measured data points for more than 150centimeter thaw depth is very much limited. Nevertheless, the results could be interpreted as Equation 4.45 could be used in Minnesota without calibration. As mentioned before, there are no data available for clayey soils in Michigan. However, since all of the developed frost and thaw depths model indicate that a statistical model developed based on the data in one state can be used in another state without calibration, data from Minnesota were used to develop a statistical model for clayey soils. As mentioned in Chapter 3, data from 5 stations were available in ten-year period. Therefore, the data from four of those stations (Ada, Marshal, Starbucks, and Rochester) were used for developing the model and data from Gatzke were used for validation. The statistical analysis is resulted in Equation 4.46 𝑇 = 8 ∗ 10−6 ∗ (𝐶𝑇𝐷𝐷)3 − 0.0051 ∗ (𝐶𝑇𝐷𝐷)2 + 1.15 ∗ 𝐶𝑇𝐷𝐷 Equation 4.46 Where all other parameters are the same as before. The results of Equation 4.46 are depicted in Figure 4.55. The calculated coefficient of determination (R2) is 0.85. Figure 4.56 shows the calculated thaw depths using Equation 4.46 versus the measured ones. The solid straight line in the figure is the locus of equality between the measured and calculated data. As can be seen, the majority of the calculated thaw depth data are within a few centimeters from the measured values. 166 Clayey soil 180 T = 8E-06CTDD3 - 0.0051CTDD2 + 1.15CTDD R² = 0.8566 n= 288 160 Thaw depth 140 120 100 80 60 40 20 0 0 50 100 150 200 250 300 o Cumulative thawing degree day ( C-day) 350 400 Figure 4.55 Thaw depths versus cumulative thawing degree-day for clayey soils in the state of Minnesota. Line of equality between the measured and calculated data Clayey soil-Minnesota Calculated thaw depth (cm) 160 140 n= 288 120 100 80 60 40 20 0 0 20 40 60 80 100 120 Measured thaw depth (cm) 140 160 Figure 4.56 Measured thaw depths in clayey soil in the state of Minnesota (Ada, Marshal, Starbucks, and Rochester) versus the thaw depth values calculated using Equation 4.46. 167 As mentioned before, to evaluate the accuracy of Equation 4.46, this equation was used to predict the measured thaw depths in Gatzke stations in the State of Minnesota. Figure 4.57 depicts the results. It should be noted that the number of measured data points in this figure is 61. Examinations of Figure 4.57 indicates that the prediction of thaw depth data in Gatzke station using Equation 4.46 is even more accurate than in the other four stations. In fact, for the calculated coefficient of determination (R2) is 0.86. Line of equality between the measured and calculated data Clayey soil-Minnesota Calculated thaw depth (cm) 160 140 n= 60 120 100 80 60 40 20 0 0 20 40 60 80 100 120 Measured thaw depth (cm) 140 160 Figure 4.57 Measured thaw depths in clayey soil in Gatzke station versus the thaw depth values calculated using Equation 4.46. In order to further evaluate the performance of Equation 4.46, performance metrics are shown in Table 4.18 for both databases. As expected, performance of Equation 4.46 indicates that the equation predicts thaw depth for both database telatively accurately. 168 Spring Load Restrictions As stated before, the accuracy of the SLR posting and removing times is critical in avoiding pavement damage. Even few days could lead to substantial damages. Most studies so far, used thawing index critical threshold to post the SLR and another critical threshold or a fixed period of 8-weeks for removing the SLR (See Chapter 2, Section 2.3.2). However, since the severity of winter varies substantially from one year to the next, using thaw depth threshold instead of thawing index threshold or a fixed period for posting and removing SLR could lead to more accurate SLR policy. Table 4.18 Performance metrics of Equation 4.46 for thaw depth estimations in the State of Minnesota Location Minnesota (Ada, Marshal, Starbucks, and Rochester) Minnesota (Gatzke) Standard Error of the Estimate (SEE) Mean Absolute Error (MAE) Mean Absolute Percentage Error (MAPE) 5.4 4.0 0.25 2.7 4.7 0.22 The Minnesota Department of Transportation (MnDOT) study showed that pavement is in the most critical state when thaw is between 15 to 30 centimeters (Ovik et.al, 2000). Their study on the rate of strength recovery and back calculation of resilient moduli also showed that at two weeks past the end of thaw the pavement strength recovers between 50 to 100 percent depends on the soil type and the increase in the fine contents in the base material led to longer recovery period (Ovik et.al, 2000). Therefore, in this study, it is recommended to post the SLR when thaw depth is at 15 centimeters and remove the SLR two weeks after the thaw completion in sandy soils and three 169 weeks after the thaw completion in clayey soils. It should be noted that these recommendations are based on literature and there was no data available to verify them. As stated before, the accuracy of the SLR posting and removing times is critical in avoiding pavement damage. On the other hand, early posting and late removal of SLR increases the trucking industry cost. Therefore, it is important to evaluate the accuracy of the thaw depth model relative to time. In other word, it is of interest to investigate the time gap between the calculated thaw depth and the measured ones. Therefore, since data in Michigan were used to develop the thaw depth model for sandy soil, data from Minnesota were used to investigate the time gap. Also, for the same reason data from Gtzake station were used to evaluate the accuracy of the thaw depth model relative to time in clayey soils. Table 4.19 shows the average time gap between the calculated and measured thaw depths in these stations. Table 4.19 Time gap between the measured and calculated thaw depth for sandy and clayey soil in the state of Minnesota Location Minnesota Sandy Soils Minnesota Clayey Soils (Gatzke) Average Absolute Time gap for 6” thaw depth (day) ±2 Average Absolute Time gap for thaw completion (day) ±3.8 Average Absolute Time gap for all thaw depth (day) ±3 ±3.6 ±8.5 ±8 As can be seen in Table 4.19, while the thaw depth model (Equation 4.46) predictions for clayey soils were more accurate than thaw depth model (Equation 4.45) for sandy soils, Equation 4.45 predictions lead to much lower time gap between the measured and calculated thaw depth. However, both models predict the beginning of the thaw with less than 4-day time gap. 170 Therefore, the accuracy of Equation 4.46 in predicting the removal of SLR should be further evaluated as more data become available. 171 CHAPTER 5 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 5.1 Summary 5.1.1 Frost Depth Modeling Frost depth is an important factor that affects the design of all infrastructures including pavements, retaining structures, building and bridge foundations and/or utility lines. Hence, accurate estimation of frost depth data plays significant roles in estimating heave pressure and heave of various structures including pavements. In this study, existing frost depth prediction models were scrutinized. These include the finite element based model (UNSAT-H) and various semi-empirical models (Stefan, Modified Berggren and Chisholm and Phang). Unfortunately, none of the model predicted the measured frost depth with and reasonable degree of accuracy. Further, some of the required input data to the models are not available and/or expensive to obtain by State Highway Agencies. Therefore, during the study, new statistical models were developed based on available and easily measured data to predict frost depths. First, the measured frost depth data in the State of Michigan were divided into two groups according to the soil types; clayey and sandy soils. For each soil type, a statistical model was developed relating the frost depth to the calculated cumulative freezing degree day (CFDD). The two statistical models were then validated using the measured frost depth data in the State of Minnesota. Both models produced reasonable estimates of the measured frost depth data. The two statistical models were then combined into one statistical model based on the average laboratory measured thermal conductivity of saturated clayey and sandy soil samples obtained from MDOT. The accuracy of the combined statistical model was then assessed using the frost 172 depth data measured in the states of Michigan and Minnesota. The calculated frost depth data were reasonable and closely represent the measured frost depth data in both states. 5.2 Frost Heave Model As soils freeze, water migrates through the soil voids below the freezing zone toward the freezing front, coats existing ice lenses causing them to grow and producing excessive frost heave. Frost have can be mitigated by removing and replacing the frost susceptible soil by drainable materials, stopping water flow by intercepting its path using drainage lines, cutting off the source of water, and/or reducing the frost depth by installing insulation. In this study, frost mitigation semi-empirical model was developed based on heat balance in the soil layer and on the newly developed statistical frost depth model. The frost mitigation model estimates the required insulation thickness to reduce or eliminate frost depth. In addition, the Gilpin’s mechanistic- empirical model was used to predict frost heave. The model yielded unreasonable results and did not simulate the frost depth data measured by MDOT. Consequently, the model was modified by replacing the heat balance equation for calculating the frost depths by the newly developed statistical frost depth model. The modified Gilpin’s model yielded relatively accurate results that represent the frost heave data measured at 5 sites in Oakland County, Michigan. Lastly, results of the frost heave model was used and heave pressure models were developed to estimate heave pressure for four soil types. 5.2.1 Thaw Depth Model During the thawing period, pavement structures, in general, are in critical conditions. The melted water at the top of the frozen area saturates the top part of the upper pavement layer causing substantial decreases in their stiffness and their load bearing capacity leading to premature and localized failure. In general, spring thaw damages are observed along city streets, 173 county roads and some state roads. These roads are designed and constructed on relatively low permeability roadbed soil and have relatively thin base and/or subbase layers. Stated differently, because of limited budget, the roadbed soils of these roads are not properly protected from frost heave. To decrease the spring thaw damage, Spring Load Restrictions (SLR) signs are usually placed along the roads. The problem stems from the fact that the time for placing and removing the SLR significantly increases the cost of pavement preservation and user costs especially the trucking industry. Hence, accurate prediction of the dates of posting and removing the SLR becomes very important. In this study, the Nixon and McRoberts thaw depth prediction model was evaluated relative to the thaw depth data measured by MDOT. The results were not satisfactory. Therefore, the cumulative thaw degree-day (CTDD) were calculated using the pavement surface and air temperatures data collected by MDOT and accurate thaw depth prediction models were developed for clayey and for sandy soils. The two models were then verified using the calculated CTDD values and thaw depth data collected by MnDOT. Finally, based on the results of the thaw depth prediction models a new policy for posting and removing SLR signs was proposed. 5.3 Conclusions Based on the results of the analyses, the following conclusions were drawn: 1. The UNSAT-H numerical finite element model for frost depth prediction requires various meteorological and soil properties data that were not available and/or expensive to collect. 2. Existing mechanistic and empirical-mechanistic models for predicting frost depths do not accurately predict the measured frost depth data in the States of Michigan. The models assume that volumetric heat capacity and water movement can be neglected. Existing statistical models are not reliable and require substantial calibration for each region. 174 3. The newly developed two statistical models for clayey and sandy soils predicted the measured frost depth data in Minnesota relatively accurately. 4. The single statistical model developed based on the average thermal conductivity of saturated clayey and sandy soils produced accurate results for both soils in the states of Michigan and Minnesota. 5. The Gilpin frost heave model was modified. The modified model yielded frost heave data that are representative to the measured data under the shoulder and under the pavement in Michigan. 6. Heave pressure model was developed based on the result of the frost heave model. However, since heave pressure data were not available, the accuracy of the model was not evaluated. 7. A new model for estimating the cumulative thawing degree day (CTDD) was developed. Based on the measured pavement surface and air temperatures in Michigan. 8. The existing mechanistic model for predicting thaw depths did not produce acceptable results. 9. Two statistical models were developed; one for clayey and one for sandy soils using the measured thaw depth data in Michigan and Minnesota. The two models predicted the measured frost depth data in Minnesota relatively accurately. 10. The two models for predicting thaw depths were used to develop a new SLR policy, the results showed that for clayey soil the model leads to average time gap of ±3 days between the measured and calculated thaw depth and for sandy soils the models leads to average time gap of ±8 days. 175 5.4 Recommendation Based on the results of this study and the enumerated conclusions, the following recommendations were made: 1. Undisturbed soil samples be collected from various soil types in Michigan and their thermal conductivity be measured. The resulting data be used to improve the accuracy of the statistical frost and thaw depth models developed in this study. 2. The developed frost model be implemented to calculate frost depth data in those areas where no temperature sensors are installed with depth. 3. The Mechanistic-Empirical Pavement Design Guide (MEPDG) or FWD data be used to estimate the range and variability of the resilient modulus of the roadbed soil and determine the sensitivity of the newly proposed SLR policy. 176 APPENDICES 177 APPENDIX A Frost and Thaw Depth Analysis 178 This appendix houses the details of frost depth calculations for the winter of 2010-2011 using Stefani equation and Modified Berggren equation for all RWIS stations in Michigan. Also provides the details of thaw depth calculations using Nixon and McRoberts equation for all of the stations. It should be noted that LP stands for Lower Peninsula and UP stands for Upper Peninsula. 179 Table A.1 Frost depth calculation using Stefan equation, Benzonia, LP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 50288 11 12 Sand 20 0.09 1.8 2.34 57924 433 444 3 γd = unit weight (Kn/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.2 Frost depth calculation using Modified Berggren equation, Benzonia, LP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 894 1.49 14788 0 894 0.00 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 50288 35499 1043 0.17 0.56 0.4 91.9 Sand 20 0.09 0.6 1565 2.42 57924 48351 1341 0.16 0.58 1.6 418.8 3 o d = layer depth (m); C= volumetric heat capacity(kJ/ m); k= thermal conductivity (W/(m. C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 82 457 0 82 457 α=v0/vs 2.9 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 9.1 vs= n(CFI)/t 3.1 FI 445 Table A.3 Thaw depth calculation using Nixon and McRoberts equation, Benzonia, LP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 o FI=freezing Index ( C-day) k 1.44 2.51 180 L 0 35201 TI 0 6 ∑TI 0 7 Table A.4 Frost depth calculation using Stefan equation, Cadillac, LP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 50288 11 11 Sand 20 0.09 2.3 2.34 57924 726 737 3 γd = unit weight (kN/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.5 Frost depth calculation using Modified Berggren equation, Cadillac, LP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 50288 35499 1229 0.17 0.56 0.4 66.3 Sand 19 0.09 1.1 1565 2.42 57924 51591 1453 0.16 0.58 1.6 746.3 3 o d = layer depth (m); C= volumetric heat capacity(kJ/ m); k= thermal conductivity (W/(m. C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 59 722 0 59 722 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 6.3 vs= n(CFI)/t 4.7 FI 729 Table A.6 Thaw depth calculation using Nixon and McRoberts equation, Cadillac, LP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 TI=Thawing Index (oC-day) 181 k 1.44 2.51 L 0 35201 TI 0 6 ∑TI 0 6 α=v0/vs 1.35 Table A.7 Frost depth calculation using Stefan equation, Grayling, LP Material γd w d k L FI ∑FI HMA 22 0 0.1 55.64 0 0 0 Gravel 20 0.08 0.3 189.00 668 10 11 Sand 19 0.09 2.1 217.73 623 633 644 3 γd = unit weight (kN/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.8 Frost depth calculation using Modified Berggren equation, Grayling, LP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 14805 0 1490 0.00 0.00 0.3 7 Gravel 20 0.08 0.3 1118 2.60 50288 39851 1227 0.17 0.56 0.4 49.8 Sand 19 0.09 1.2 1565 2.42 57931 53114 1475 0.16 0.58 1.7 677.7 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 50 653 0 50 653 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 4.5 vs= n(CFI)/t 5.2 FI 637 Table A.9 Thaw depth calculation using Nixon and McRoberts equation, Grayling, LP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 Sand 19 1.08 0.1 TI=Thawing Index (oC-day) 182 k 1.44 2.51 2.34 L 0 35201 40565 TI 0 6 20 ∑TI 0 6 29 α=v0/vs 0.86 Table A.10 Frost depth calculation using Stefan equation, Houghton Lake, LP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 50288 11 12 Loose Sand 20 0.09 2.1 2.34 57924 618 629 3 γd = unit weight (kN/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.11 Frost depth calculation using Modified Berggren equation, Houghton Lake, LP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 50288 39858 1229 0.17 0.56 0.4 50.0 Sand 20 0.09 1.2 1565 2.42 57924 53007 1453 0.16 0.58 1.7 634.4 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 52 616 0 52 616 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 4.5 vs= n(CFI)/t 5.2 FI 629 Table A.12 Thaw depth calculation using Nixon and McRoberts equation, Houghton Lake, LP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 Loose Sand 20 1.08 0.1 TI=Thawing Index (oC-day) 183 k 1.44 2.51 L 0 35201 TI 0 18 ∑TI 0 19 2.34 40565 35 53 α=v0/vs 0.84 Table A.13 Frost depth calculation using Stefan equation, Ludington, LP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 50288 11 12 Sand with Clay 18 0.09 1.2 2.34 92529 344 356 3 γd = unit weight (kN/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.14 Frost depth calculation using Modified Berggren equation, Ludington, LP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 50288 39858 1229 0.17 0.56 0.4 58.1 Sand 18 0.09 0.5 1565 2.42 92529 68950 1416 0.16 0.58 0.8 357.5 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 62 369 0 62 369 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 5.1 184 vs= n(CFI)/t 3.1 FI 361 α=v0/vs 1.63 Table A.15 Frost depth calculation using Stefan equation, Reed City, LP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 52299 11 11 Sand 20 0.09 1.9 2.34 60345 526 536 3 γd = unit weight (kN/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.16 Frost depth calculation using Modified Berggren equation, Reed City, LP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 52299 41273 1229 0.17 0.56 0.4 53.1 Sand 20 0.09 1.0 1565 2.42 60345 54497 1453 0.16 0.58 1.4 550.6 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 56 545 0 56 545 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 4.7 vs= n(CFI)/t 4.6 FI 533 Table A.17 Thaw depth calculation using Nixon and McRoberts equation, Reed City, LP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 TI=Thawing Index (oC-day) 185 k 1.44 2.51 L 0 35201 TI 0 16 ∑TI 0 17 α=v0/vs 1.03 Table A.18 Frost depth calculation using Stefan equation, Waters, LP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 50288 11 11 Sand 20 0.09 0.6 2.51 50288 42 52 Loose Sand with Gravel 19 1.09 2.1 2.34 57924 603 655 γd = unit weight (kN/m3); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.19 Frost depth calculation using Modified Berggren equation, Waters, LP Material γd w d C k L Ĺ= ∑Ld/∑d 0 39858 45967 Ĉ= ∑Cd/∑d 1490 1229 1155 μ= Ĉ / Ĺ *vs 0.00 0.17 0.16 λ R ∑R+R/2 FI ∑FI 0 51 233 0 51 233 648 648 α=v0/vs 0.77 HMA 22 0 0 1490 1.49 0 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 50288 0.56 0.4 48.8 Sand 20 0.09 0.6 1118 2.60 50288 0.58 0.8 205.0 Loose Sand with 19 1.09 0.7 1565 2.42 57924 50660 1304 0.16 0.58 1.0 466.9 Gravel d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 4.1 5.3 vs= n(CFI)/t FI 656 Table A.20 Thaw depth calculation using Nixon and McRoberts equation, Waters, LP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 Sand 20 1.08 0.1 TI=Thawing Index (oC-day) 186 k 1.44 2.51 2.51 L 0 36617 35201 TI 0 14 6 ∑TI 0 14 19 Table A.21 Frost depth calculation using Stefan equation, Williamsburg, LP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 50288 10 11 Sand 20 0.09 0.3 2.34 57924 10 21 Silty Clay 18 1.09 1.1 1.83 101767 388 409 γd = unit weight (kN/m3); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.22 Frost depth calculation using Modified Berggren equation, Williamsburg, LP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.4 1118 2.60 50288 35499 1229 0.17 0.56 0.4 55.6 Sand 20 0.09 0.7 1453 2.42 57924 44216 1341 0.16 0.58 0.4 100.6 Silty Clay 18 1.09 1.1 1565 1.90 101767 64890 1527 0.16 0.58 0.7 286.3 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 49 139 393 0 49 139 393 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 5.1 187 vs= n(CFI)/t 3.4 FI 408 α=v0/vs 1.45 Table A.23 Frost depth calculation using Stefan equation, Au Train, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 52299 11 11 Sand 20 0.09 0.5 2.34 60345 34 46 Loose Sand 19 1.09 2.7 2.34 57924 1027 1072 γd = unit weight (kN/m3); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.24 Frost depth calculation using Modified Berggren equation, Au Train, UP Material γd w d C k L Ĺ= ∑Ld/∑d 0 36915 49356 Ĉ= ∑Cd/∑d 1490 1229 1416 μ= Ĉ / Ĺ *vs 0.00 0.17 0.16 λ R ∑R+R/2 HMA 22 0 0 1490 1.49 0 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 52299 0.56 0.4 49.4 Sand 20 0.09 0.5 1565 2.42 60345 0.58 0.7 186.3 Loose 19 1.09 1.2 1565 2.42 57924 54236 1490 0.16 0.58 1.8 966.3 Sand d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; FI ∑FI 0 44 209 0 44 209 1068 1068 α=v0/vs 0.44 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 3.8 vs= n(CFI)/t 8.8 FI 1072 Table A.25 Thaw depth calculation using Nixon and McRoberts equation, Au Train, UP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 Sand 19 1.08 0.1 TI=Thawing Index (oC-day) 188 k 1.44 2.51 2.34 L 0 36617 42242 TI 0 15 3 ∑TI 0 14 18 Table A.26 Frost depth calculation using Stefan equation, Brevort, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 52299 11 11 Sand 20 0.09 2.2 2.34 60345 707 718 3 γd = unit weight (kN/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.27 Frost depth calculation using Modified Berggren equation, Brevort, UP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 52299 36915 1229 0.17 0.56 0.4 41.3 Sand 20 0.09 1.3 1565 2.42 60345 54534 1490 0.16 0.58 1.9 761.9 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 37 714 0 37 714 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 3.7 vs= n(CFI)/t 5.3 FI 711 Table A.28 Thaw depth calculation using Nixon and McRoberts equation, Brevort, UP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 Sand 20 1.08 0.1 TI=Thawing Index (oC-day) 189 k 1.44 2.51 2.34 L 0 35201 40565 TI 0 18 3 ∑TI 0 18 21 α=v0/vs 0.7 Table A.29Frost depth calculation using Stefan equation, Cooks, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 52299 11 11 Sand 20 0.09 1.9 2.34 60345 525 536 3 γd = unit weight (kN/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.30 Frost depth calculation using Modified Berggren equation, Cooks, UP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 50288 35499 1229 0.17 0.56 0.4 46.9 Sand 20 0.09 1.1 1565 2.42 57924 51480 1416 0.16 0.58 1.6 545.6 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 42 526 0 42 526 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 3.7 190 vs= n(CFI)/t 3.7 FI 533 α=v0/vs 1.03 Table A.31 Frost depth calculation using Stefan equation, Engadine, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 52299 11 11 Sand 20 0.09 0.6 2.34 60345 58 69 Silty Clay 18 1.09 1.4 1.83 101767 551 620 γd = unit weight (kN/m3); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.32 Frost depth calculation using Modified Berggren equation, Engadine, UP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.4 1118 2.60 52299 36915 1229 0.17 0.56 0.4 43.1 Sand 20 0.09 0.7 1453 2.42 60345 50623 1341 0.16 0.58 0.9 228.8 Silty Clay 18 1.09 1.1 1565 1.90 101767 70105 1490 0.16 0.58 1.0 415.0 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 38 242 611 0 38 242 611 α=v0/vs 0.66 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 3.8 191 vs= n(CFI)/t 5.7 FI 617 Table A.33 Frost depth calculation using Stefan equation, Golden Lake, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 50288 10 11 Sand 20 0.09 2.7 2.34 57924 989 999 3 γd = unit weight (kN/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.34 Frost depth calculation using Modified Berggren equation, Golden Lake, UP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 50288 35499 1229 0.17 0.56 0.4 43.8 Sand 20 0.09 1.8 1453 2.42 57924 53454 1490 0.16 0.58 2.5 1068.1 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 39 988 0 39 988 α=v0/vs 0.39 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 2.9 vs= n(CFI)/t 7.4 FI 1002 Table A.35 Thaw depth calculation using Nixon and McRoberts equation, Golden Lake, UP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 Sand 20 1.08 0.1 TI=Thawing Index (oC-day) 192 k 1.44 2.51 2.34 L 0 36617 42242 TI 0 19 3 ∑TI 0 19 22 Table A.36 Frost depth calculation using Stefan equation, Harvey, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 55.64 0 0 0 Gravel 20 0.08 0.3 196.56 668 11 12 Sand with Gravel and Silt 20 0.09 0.3 226.80 623 18 29 Sand 20 1.09 2.0 189.00 668 453 483 γd = unit weight (kN/m3); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oCday) Table A.37 Frost depth calculation using Modified Berggren equation, Harvey, UP Material γd w d C k L Ĺ= ∑Ld/∑d 0 36915 Ĉ= ∑Cd/∑d 1490 1229 μ= Ĉ / Ĺ *vs 0.00 0.17 λ R ∑R+R/2 HMA 22 0 0 1490 1.49 0 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 52299 0.56 0.4 45.6 Sand with Silt 19 0.09 0.5 1453 2.42 57924 48053 1416 0.16 0.58 0.7 166.3 Sand 20 1.09 0.6 1118 2.60 50288 48984 1304 0.16 0.58 0.9 320.6 3 o d = layer depth (m); C= volumetric heat capacity(kJ/ m); k= thermal conductivity (W/(m. C); μ= fusion parameter; FI ∑FI 0 41 0 41 188 473 188 473 α=v0/vs 0.74 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 2.9 4.0 vs= n(CFI)/t FI 480 Table A.38 Thaw depth calculation using Nixon and McRoberts equation, Harvey, UP Material γd w HMA 22 0 Gravel 20 0.08 Sand with Silt 19 1.08 o TI=Thawing Index ( C-day) d 0.1 0.3 0.1 193 k 1.44 2.51 2.34 L 0 40565 35201 TI 0 7 3 ∑TI 0 7 10 Table A.39 Frost depth calculation using Stefan equation, Michigamme, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 52299 11 11 Sand 20 0.09 0.6 2.34 60345 54 64 Clayey Sand 19 1.09 2.3 1.83 57924 904 969 γd = unit weight (kN/m3); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.40 Frost depth calculation using Modified Berggren equation, Michigamme, UP Material γd w d C k L Ĺ= ∑Ld/∑d 0 36915 50623 Ĉ= ∑Cd/∑d 1490 1229 1416 μ= Ĉ / Ĺ *vs 0.00 0.17 0.16 λ R ∑R+R/2 HMA 22 0 0 1490 1.49 0 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 52299 0.56 0.4 45.6 Sand 20 0.09 0.6 1453 2.42 60345 0.58 0.9 240.0 Clayey 19 1.09 1.3 1118 1.90 57924 54646 1267 0.16 0.58 2.4 783.8 Sand d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; FI ∑FI 0 41 254 0 41 254 951 951 α=v0/vs 0.5 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 3.2 194 vs= n(CFI)/t 6.5 FI 964 Table A.41 Frost depth calculation using Stefan equation, Seney, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 52299 11 11 Sand 20 0.09 2.1 2.34 60345 624 619 3 γd = unit weight (kN/m ); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.42 Frost depth calculation using Modified Berggren equation, Seney, UP Ĺ= Ĉ= μ= λ R ∑R+R/2 ∑Ld/∑d ∑Cd/∑d Ĉ / Ĺ *vs HMA 22 0 0 1490 1.49 0 0 1490 0.00 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 52299 36915 1229 0.17 0.56 0.4 51.9 Sand 20 0.09 1.2 1565 2.42 60345 53975 1490 0.16 0.58 1.7 658.1 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; Material γd w d C k L FI ∑FI 0 46 631 0 46 631 α=v0/vs 0.78 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 4.1 vs= n(CFI)/t 5.2 FI 618 Table A.43 Thaw depth calculation using Nixon and McRoberts equation, Seney, UP Material γd w d HMA 22 0 0.1 Gravel 20 0.08 0.3 Sand 20 1.08 0.2 TI=Thawing Index (oC-day) 195 k 1.44 2.51 2.34 L 0 36617 42242 TI 0 9 11 ∑TI 0 9 20 Table A.44 Frost depth calculation using Stefan equation, St. Ignace, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 52299 11 11 Sand 20 0.09 0.6 2.34 60345 54 64 Silty Clayey Sand 18 1.09 1.3 1.83 101767 542 607 γd = unit weight (kN/m3); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.45 Frost depth calculation using Modified Berggren equation, St. Ignace, UP Material γd w d C k L Ĺ= ∑Ld/∑d 0 41273 52448 Ĉ= ∑Cd/∑d 1490 1229 1341 μ= Ĉ / Ĺ *vs 0.00 0.17 0.16 λ R ∑R+R/2 HMA 22 0 0 1490 1.49 0 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 52299 0.56 0.4 41.3 Sand 20 0.09 0.6 1453 2.42 60345 0.58 0.9 213.1 Silty Clayey 10176 18 1.09 0.6 1788 1.90 71222 1527 0.16 0.58 1.2 406.3 Sand 7 d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; FI ∑FI 0 37 226 0 37 226 587 587 α=v0/vs 0.4 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 2.2 196 vs= n(CFI)/t 5.2 FI 600 Table A.46 Frost depth calculation using Stefan equation, Twin Lakes, UP Material γd w d k L FI ∑FI HMA 22 0 0.1 1.44 0 0 0 Gravel 20 0.08 0.3 2.51 52299 11 11 Sand 20 0.09 0.6 2.34 60345 54 64 Silty Clayey Sand 18 1.09 1.6 1.83 101767 796 861 γd = unit weight (kN/m3); w = water content (%); d = layer depth (m); k= thermal conductivity (W/(m.o C); L= latent heat of fusion (kJ/3m); FI=freezing Index (oC-day) Table A.47 Frost depth calculation using Modified Berggren equation, Twin Lakes, UP Material γd w d C k L Ĺ= ∑Ld/∑d 0 36915 50623 Ĉ= ∑Cd/∑d 1490 1229 1341 μ= Ĉ / Ĺ *vs 0.00 0.17 0.16 λ R ∑R+R/2 HMA 22 0 0 1490 1.49 0 0.00 0.3 0 Gravel 20 0.08 0.3 1118 2.60 52299 0.56 0.4 41.3 Sand 20 0.09 0.6 1453 2.42 60345 0.58 0.9 218.1 Silty Clayey 18 1.09 0.9 1788 1.90 101767 74128 1565 0.16 0.58 1.7 688.8 Sand d = layer depth (m); C= volumetric heat capacity(kJ/3m); k= thermal conductivity (W/(m.o C); μ= fusion parameter; FI ∑FI 0 37 231 0 37 231 843 843 α=v0/vs 0.5 λ= correction factor; R= thermal diffusivity (m2.oC/W); FI=freezing Index (oC-day) v0= Annual Temp average-0 3.1 6.3 vs= n(CFI)/t FI 881 Table A.48 Thaw depth calculation using Nixon and McRoberts equation, Twin Lakes, UP Material γd w HMA 22 0 Gravel 20 0.08 Sand 20 1.08 o TI=Thawing Index ( C-day) d 0.1 0.3 0.1 197 k 1.44 2.51 2.34 L 0 35201 42242 TI 0 9 3 ∑TI 0 9 12 APPENDIX B Frost Heave Stations Profile 198 This appendix houses the details pavement profile of frost heave stations in Michigan. In the figures both measured frost heave and frost depth were shown for each station (Novak, 1968). 199 1963 1964 Heave Rate at Pavement and Shoulder Edge, cm 1.3 2.5 1.3 2.5 Frozen Zone Below Pavement and Shoulder, cm 0 13 25 38 51 64 76 89 PAVEMENT UNDERCUT 30 CM FOR FROST PROTECTION 102 ' Figure B.1 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 528+88. 200 1963 1964 Heave Rate at Pavement and Shoulder Edge, cm 1.3 2.5 1.3 2.5 Frozen Zone Below Pavement and Shoulder, cm 0 13 25 38 51 64 76 89 102 Figure B.2 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 652+00. 201 1963 1964 Heave Rate at Pavement and Shoulder Edge, cm 1.3 2.5 0 1.3 2.5 Frozen Zone Below Pavement and Shoulder, cm 0 13 25 38 51 64 76 89 102 Figure B.3 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 719+00. 202 1964 Heave Rate at Pavement and Shoulder Edge, cm 1963 1.3 2.5 1.3 2.5 Frozen Zone Below Pavement and Shoulder, cm 0 13 25 38 51 64 76 89 102 Figure B.4 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 724+00. 203 1963 1964 1.3 Heave Rate at Pavement and Shoulder Edge, cm 2.5 0 1.3 2.5 Frozen Zone Below Pavement and Shoulder, cm 0 13 25 38 51 64 76 89 102 Figure B.5 Frost depth and corresponding frost heave of shoulder and pavement, Sta. 474+00. 204 REFERENCES 205 REFERENCES 1. American Concrete Pavement Association, “Frost-Susceptible Soils”, 2008. 2. Aldrich Jr, P. Harl, and H. M. Paynter, “Analytical Studies of Freezing and Thawing of Soils”, Technical Report No. 42, Arctic Construction and Frost Effects Laboratory, New England Division U.S. Army Corps of Engineers, Boston, MA, U.S. Arctic Construction and Frost Effects Lab Boston Ma, 1953. 3. Andersland, O. B., and D. M Anderson, “Geotechnical Engineering for Cold Regions”, McGraw-Hill Inc, 1978. 4. ASTM C578-14a., “Standard Specification for Rigid, Cellular Polystyrene Thermal Insulation”, ASTM International, West Conshohocken, PA, 2014, www.astm.org. 5. Berg, R. L., “Calculating Maximum Frost Depths At Mn/Road: Winters 1993-94, 1994-95 And 1995-96”, No. MN/RC-97/21, 1997. 6. Bianchini, A., and R.G. Carlos., “Pavement-Transportation Computer Assisted Structural Engineering (PCASE) Implementation of the Modified Berggren (ModBerg) Equation for Computing the Frost Penetration Depth within Pavement Structures”, Publication No. ERDC/GSL-TR-12-15. Engineer Research and Development Center Vicksburg Ms Geotechnical and Structures Lab, 2012. 7. Black, P. B., “Applications of the Clapeyron Equation to Water and Ice in Porous Media”, No. CRREL-95-6. Cold Regions Research and Engineering Lab Hanover, 1995. 8. Boselly III, S. E., G. S. Doore, J. E. Thornes, C. Ulberg, and D. Ernst, “Road Weather Information Systems Volume 1”, Research report. Strategic Highway Research Program Publication-SHRP-H-350, National Research Council, Washington, DC, 1993 9. Boyd, D. W., “Normal Freezing and Thawing Degree-Days from Normal Monthly Temperatures” Canadian Geotechnical Journal, Volume 13.2, pp.176-180, 1976. 10. Bradley, H., M. A. Ahammed, S. Hilderman, and S. Kass, “Responding to Climate Change with Rational Approaches for Managing Seasonal Weight Programs in Manitoba.” Cold Regions Engineering, pp. 391-401, 2012. 11. Chapin, J., B. H. Kjartanson, and J. Pernia, “Comparison of Two Methods for Applying Spring Load Restrictions on Low Volume Roads”, Canadian Journal of Civil Engineering, Vol. 39, Issue 6, pp.599-609, 2012. 12. Chisholm, R. A., and W. A. Phang, “Measurement and Prediction of Frost Penetration in Highway”, Transportation Research Record: Journal of Transportation Research Board, No 205 918, Transportation Research Board of the National Academies, Washington D.C., pp. 1-10, 1983. 13. Crandell, J. H., "Below-Ground Performance of Rigid Polystyrene Foam Insulation: Review of Effective Thermal Resistivity Values Used in ASCE Standard 32-01—Design and Construction of Frost-Protected Shallow Foundations”, Journal of Cold Regions Engineering, Volume 24.2, pp. 35-53, 2009. 14. Dash, J. G., A. W. Rempel, and J. S. Wettlaufer, “The Physics of Premelted Ice and Its Geophysical Consequences”, Reviews of Modern Physics, Vol. 78.3, pp. 695-741, 2006. 15. Departments of The Army and The Air Force, “Artic and Subarctic Construction Calculation Methods for Determination of Depths of Freeze and Thaw in Soil”, Technical Manual, No. 5852-6, Air Force Regulation, AFR88-19, Volume 6, January, 25, 1988. 16. Devices, Decagon, “KD2 Pro Thermal Properties Analyzer Operator's Manual, Version 10.0”, 2008. 17. Diefenderfer, B. K., Imad L. Al-Qadi, and Stacey D. Diefenderfer, “Model to Predict Pavement Temperature Profile: Development and Validation”, Journal of Transportation Engineering, Vol. 132.2, 2006, pp. 162-167. 18. Everett, D. H., “The Thermodynamics of Frost Damage to Porous Solids”, Transactions of the Faraday Society, Volume 57, pp. 1541-1551, 1961. 19. Edgar, T. V., “Instrumentation and Analysis of Frost Heave Mitigation on WY-70, Encampment, WY”.No. FHWA-WY-14/01F, 2014. 20. Fowler, A. C., and W. B. Krantz, “A Generalized Secondary Frost Heave Model”, Siam Journal on Applied Mathematics, Volume 54.6, pp. 1650-1675, 1994. 21. Fayer, M.l J., “UNSAT-H Version 3.0: Unsaturated Soil Water and Heat Flow Model. Theory, User Manual, and Examples”, Pacific Northwest National Laboratory 13249, 2000. 22. Gilpin, R. R., “A Model of the 'Liquid-Like' Layer between Ice and a Substrate With Applications to Wire Regelation and Particle Migration”, Journal of Colloid Interface Science, Volume 68, pp.235-251, 1979. 23. Gilpin, R. R., “Theatrical Studies of Particle Engulfment”, Journal of Colloid Interface Science, Volume 74, pp.44-63, 1980 24. Gilpin, R. R., “A Model for the Prediction of Ice Lensing and Frost Heave in Soils”, Water Resources Research, Volume 16.5, pp. 918-930, 1980. 25. Gilpin, R. R., “A Frost Heave Interface Condition for Use in Numerical Modeling”, Proc. 4th Canadian Permafrost Conference, Calgary, Canada, pp. 459-465, 1982. 206 26. Han, S. J., and D. J. Goodings, “Practical Model of Frost Heave in Clay”, Journal of Geotechnical and Geoenvironmental Engineering, Volume 132.1, pp.92-101, 2006. 27. Harlan, R. L., “Analysis of Coupled Heat‐Fluid Transport in Partially Frozen Soil”, Water Resources Research, Volume 9.5, pp. 1314-1323, 1973. 28. Hayhoe, H. N., and D. Balchin, “Field frost heave measurement and prediction during periods of seasonal frost.” Canadian Geotechnical Journal, Volume 27.3, pp. 393-397, 1990. 29. Hsieh, C. K., C. Qin, and E. E. Ryder, “Development of Computer Modeling for Prediction of Temperature Distribution inside Concrete Pavements”, Report No. FL/DOT/SMO/90-374, Florida Department of, 1989. 30. Huen, K., S. Tighe, B. Mills, and M. Perchanok, “Development of Tools for Improved Spring Load Restriction Policies in Ontario”, The Annual Conference of The Transportation Association of Canada (TAC), Charlottetown, Prince Edward Island, 2006. 31. Iqbal, M., “An Introduction to Solar Radiation”, Elsevier, 2012. 32. Janoo, V., and E. Cortez, “Pavement Evaluation in Cold Regions”, Cold Regions Engineering, pp.360-371, 2002. 33. Jiji, L. M., “Heat Conduction”, Springer, 2009. 34. Klein, J., “Finite element method for time‐dependent problems of frozen soils”, International Journal for Numerical and Analytical Methods in Geomechanics, Volume 5.3, pp. 263-283, 1981. 35. Kestler, M. A., R. L. Berg, B. C. Steinert, G. L. Hanek, M. A. Truebe, and D. N. Humphrey, “Determining When to Place and Remove Spring Load Restrictions on Low-Volume Roads: Three Low-Cost Techniques”, Transportation Research Record: Journal of the Transportation Research Board, Volume 1, pp. 219-229, 2007. 36. Kinosita, S., “Heaving force of frozen soils”, Physics of Snow and Ice: proceedings, Volume 1.2, pp. 1345-1360, 1967. 37. Konrad, J.M., and N. R. Morgenstern, “A Mechanistic Theory of Ice Lens Formation in FineGrained Soils”, Canadian Geotechnical Journal, Volume 17.4, pp.473-486, 1980. 38. Konrad, J-M., and N. R. Morgenstern, “Effects of applied pressure on freezing soils”, Canadian Geotechnical Journal, Volume.19.4, pp. 494-505, 1982. 39. Lackner, R., C. Pichler, and A. Kloiber, “Artificial ground freezing of fully saturated soil: viscoelastic behavior”, Journal of engineering mechanics, Volume 134.1, pp. 1-11, 2008. 40. Lai, Y, W. Hui, W. Ziwang, L. Songyu, and D. Xuejun, “Analytical viscoelastic solution for frost force in cold-region tunnels”, Cold regions science and technology 31, Volume 3, pp. 227-234, 2000. 207 41. Leong, P., S. Tighe, and G. Doré, “Using LTPP data to develop spring load restrictions: A pilot study”, Proceeding the 84th Annual Meeting CD-ROM, Transportation Research Board, Washington, D.C., 2005. 42. Levinson, D. M., O. Mihai, V. Voller, I Margineau, B. Smalkowski, M. Hashami, N. Li, M. Corbett, and E. Lukanen, “Cost/Benefit Study of Spring Load Restrictions”, 2005. 43. Liu, Z., X. Yu, Y. Sun, and B. Zhang, “Formulation and Characterization of Freezing Saturated Soils”, Journal of Cold Regions Engineering, Volume 27(2), pp. 94-107, 2012. 44. Loch, J. P. G., and R. D. Miller, “Tests of The Concept of Secondary Frost Heaving”, Soil Science Society of America Journal, Volume 39.6, pp. 1036-1041, 1975. 45. Marquis, B., “Mechanistic Approach to Determine Spring Load Restrictions in Main”, Technical report No. 08-1, Main Department of Transportation, Bangor, Main, 2008. 46. Minnesota Department of Transportation (MnDOT), “Guidelines for Seasonal Load Limit Starting and Ending Dates”, Policy, Safety & Strategic Initiatives Division Technical Memorandum. No. 09-09-MAT-02, 2009. 47. Michigan Department of Transportation (MDOT), “Plans of Proposed”, No. NH 0852(019), pp.17-18, 2008. 48. Michigan Department of Transportation (MDOT), “Plans of Proposed”, No. ARRA 0984(084), pp.27, 2009a. 49. Michigan Department of Transportation (MDOT), “Request For Proposals, Book 2 Project Requirements”, North and Superior Regions ITS Deployment, pp.81-110, 2009b. 50. Michigan State University Enviro Weather http://www.enviroweather.msu.edu/homeMap.php, Accessed November, 2012 (MSU-EW) 51. Miller, H., C. Cabral, M. Kestler, R. Berg, and R. Eaton, “Calibration of a Freeze-Thaw Prediction Model for Spring Load Restriction Timing in Northern New England”, Cold Regions Engineering, pp. 369-379, 2012 52. Mnhoney, J., R. Rutherford and G.Hicks, “Guidelines for Spring Highway Use Restrictions”, Report No. WA-RD-80.2, Washington State Department of Transportation, Olympia, WA, 1986. 53. National Oceanic and Atmospheric Administration (NOAA). 54. Nixon, J. and E. McRoberts, “A study of Some Factors Affecting the Thawing of Frozen Soils”, Canadian Geotechnical Journal, Volume 10, pp. 439-452, 1973. 55. Nixon, J. F., “Discrete ice lens theory for frost heave in soils”, Canadian geotechnical journal, Volume 28.6, 843-859, 1991. 208 56. Novak, E. C., “Study of Frost Action in Class AA Shoulders near Pontiac”, Research report No.671, Project 63E-27, Testing and Research Division, Department of State Highway, Michigan, 1968. 57. Ovik, J.M., J.A. Siekmeier, and D.A Van Duesen, “Improved Spring Load Restriction Guidelines Using Mechanistic Analysis”, Minnesota Department of Transportation, Report No. MN/RC-2000-18, 2000. 58. Penner, E., “Frost heaving forces in Leda clay”, Canadian Geotechnical Journal, Volume 7.1 , pp. 8-16, 1969. 59. Penner, E., and L. W. Gold., “Transfer of heaving forces by adfreezing to columns and foundation walls in frost-susceptible soils”, Canadian Geotechnical Journal, Volume 8.4, pp. 514-526, 1971. 60. Penner, E., and W. W. Irwin, “Adfreezing of Leda clay to anchored footing columns”, Canadian Geotechnical Journal. Volume 6, pp. 327-337, 1969. 61. Peppin, S. S., and R. W. Style, “The Kinetics of Ice-Lens Growth in Porous Media”, Journal of Fluid Mechanics, Vol. 692, pp. 482-498, 2012. 62. Peppin, S. S., and R. W. Style, “The Physics of Frost Heave and Ice-Lens Growth”, Vadose Zone Journal, Vol.12, 2013. 63. Road Weather Information System (RWIS), “Road Weather Information System portal, Michigan”, http://www.rwisonline.com/scanweb/swlogin.asp, Accessed November , 2012. 64. Smith, J. M., “Introduction to Chemical Engineering Thermodynamic”, McGraw-Hill, Boston, 2001. 65. Stormont, J., and S. Zhou, “Improving Pavement Sub-surface Drainage Systems by Considering Unsaturated Water Flow”, U.S. Department of Transportation Federal Highway Administration, Department of Civil engineering University of New Mexico, 2001 66. Taber, S., “The Mechanics of Frost Heaving”, The Journal of Geology, Volume 2, pp. 303317, 1930. 67. Takagi, S., “Segregation Freezing as the Cause of Suction Force For Ice Lens Formation”, Engineering Geology, Volume 13.1, pp. 93-100, 1979. 68. Tighe, S. L., B. Mills, C.T. Haas., and S. Baiz, “Using Road Weather Information Systems (RWIS) to Control Load Restrictions on Gravel and Surface-Treated Highways”, Ministry of Transportation Engineering Standards Branch Report, 2007. 69. Thompson, M. R., Dempsey, B. J., Hill, H., and J. Vogel, “Characterizing Temperature Effects for Pavement Analysis and Design”, Transportation Research Record: Journal of Transportation Research Board, No 1121, Transportation Research Board of the National Academies, Washington D.C., pp. 14-22, 1987. 209 70. Transportation Officials, “AASHTO Guide for Design of Pavement Structures” ASSHTO, Volume 1, 1993. 71. US department of transportation, “An Introduction to Standards for Road Weather Information Systems (RWIS)”, Federal Highway Administration, 2002. 72. U. S. Army Corps of Engineers, “Pavement Criteria for Seasonal Frost Conditions”, Department of the Army Technical Manual, EM 1110-3-138, 1984. 73. Vialov, S. S., V. G. Gmoshinskii, S. E. Gorodetskii, V. G. Grigorieva, Iu K. Zaretskii, N. K. Pekarskaia, and E. P. Shusherina, “The strength and creep of frozen soils and calculations for ice-soil retaining structures”, No. Translation 76. 1965. 74. Wallace, M., “How to Prevent Frost Heave”, Concrete Construction, Volume 32.4, pp. 369372, 1987. 75. Wisconsin Department of Transportation, “Using Weight Limits to Protect Local Roads”, Bulletin. No. 8, 2003 76. Yavuzturk, C., K. Ksaibati, and A. D. Chiasson, “Assessment of Temperature Fluctuations in Asphalt Pavements Due to Thermal Environmental Conditions Using a Two-Dimensional, Transient Finite-Difference Approach”, Journal of Materials in Civil Engineering, Volume.17.4, pp.465-475, 2005. 77. Yoder, El. J., and M. W. Witczak, “Principles of Pavement Design”, John Wiley & Sons, 1975. 78. Yuanming, L., Xuefu, Z., Jianzhang, X., Shujuan, Z., & Zhiqiang, L., “Nonlinear analysis for frost-heaving force of land bridges on Qing-Tibet Railway in cold regions”, Journal of Thermal Stresses, Volume 28(3), pp. 317-331, 2005. 79. Zapata, C. E., and W. N. Houston, “Calibration and Validation of the Enhanced Integrated Climatic Model for Pavement Design”, Transportation Research Board, Volume 602, 2008. 80. ARA, Inc., “Guide for mechanistic-empirical design of new and rehabilitated pavement structures”, 2004. 81. Solaimanian, Mansour, and Pablo Bolzan., “Analysis of the integrated model of climatic effects on pavements”, No. SHRP-A-637. Strategic Highway Research Program, National Research Council, 1993. 82. Dempsey, Barry J., and Richard A. Pur, “Guidelines for Design, Construction, and Evaluation of Airport Pavement Drainage”, Illinois University at Urbana Department of Civil Engineering, 1990. 83. Johanneck, Luke Anton, “A comprehensive evaluation of the effects of climate in MEPDG predictions and of MEPDG EICM model using MnROAD data”, Master Dissertation., University of Minnesota, 2011. 210 84. Liang, Robert Y., “Validation of enhanced integrated climatic model prediction over different drainable base materials”, In Transportation Research Board 85th Annual Meeting, No. 06-2529. 2006. 85. Heydinger, A., “Monitoring seasonal instrumentation and modeling climatic effects on pavement at the ohio/shrp test road”, No. FHWA/OH-2003/018, 2003. 86. Ahmed, Zubair, Ivana Marukic, Sameh Zaghloul, and Nick Vitillo, “Validation of enhanced integrated climatic model predictions with New Jersey seasonal monitoring data”, Transportation Research Record: Journal of the Transportation Research Board 1913, pp. 148-161, 2005. 87. McCartney, John S., R. Panneer Selvam, Josh King, and Ali Khosravi, “Evaluation of the enhanced integrated climatic model for the Arkansas State Highway and Transportation Department”, No. TRC-0902. 2010. 88. Khazanovich, L., Balbo, J.T., Johanneck, L., Lederle, R., Marasteanu, M., Saxena, P., Tompkins, D., Vancura, M., Watson, M., Harvey, J. and Santero, N.J., “Design and construction guidelines for thermally insulated concrete pavements”, 2013. 89. Chung, Yoonseok, and Alex Hak-Chul Shin, “Local Calibration of EICM Using Measured Temperature Gradients and Numerical Analysis”, International Journal of Pavement Research and Technology 8, No. 4, pp. 259- 266, 2015. 90. Yesiller, Nazli, Craig H. Benson, and Peter J. Bosscher, “Comparison of load restriction timings determined using FHwA guidelines and frost tubes”, Journal of cold regions engineering 10.1, pp. 6-24, 1996. 91. Miller, H., C. Cabral, M. Kestler, R. Berg, and R. Eaton, “Comparative Analyses of Methods for Posting Spring Load Restrictions”, In Proceedings of the 10th International Symposium on Cold Regions Development, Anchorage, Alaska June, pp. 2-5, 2013. 92. Berg, R. L., M. A. Kestler, R. A. Eaton, and C. C. Benda, “Estimating when to apply and remove spring load restrictions”, In Proceedings of the American Society of Civil Engineers 13th International Conference on Cold Regions Engineering. CD-ROM. Orono, ME. 2006. 93. Kestler, Maureen, Richard Berg, Heather Miller, Bryan Steinert, Robert Eaton, Gregg Larson, and John Haddock, “Keeping Springtime Low-Volume Road Damage to a Minimum: Toolkit of Practical Low-Cost Methods for Road Managers”, Transportation Research Record: Journal of the Transportation Research Board, Volume 2205, pp. 155-164, 2011. 94. Doré, Guy, “Development and Validation of the Thaw-weakening Index”, International Journal of Pavement Engineering 5, No. 4, pp.185-192, 2004. 211