Validation of a dynamic modulus predictive equation on the basis of spanish asphalt concrete mixtures

A. Mateos^a*, J.B. Soares^b 

a. Centro de Estudios del Transporte del CEDEX, (Madrid, Spain)

b. Universidade Federal do Ceará, (Fortaleza, Ceará, Brazil)

*angel.mateos@cedex.es

1. INTRODUCTIONTOP

Dynamic modulus (|E*|) is one of the most important properties of asphalt mixtures, since it determines the strain response characteristics as a function of loading rate and temperature. It is defined as the ratio of peak cyclic stress to peak cyclic strain under harmonic loading. It can be used to predict pavement structural response to traffic loads (1), and it can be regarded as a performance-related material property for mixture specifications and design. The latter is the case of the European Standard EN 13108-1, which includes stiffness as one of the requirements within the fundamental approach for asphalt concrete characterization. A similar situation exists in the United States of America, where dynamic modulus is one of the candidate tests to complement Superpave volumetric mix design (2). Consequently, it is not surprising that the great majority of the existing analytical methods for pavement design require the determination of this parameter. The most advanced generation of analytical methods, the incremental-recursive procedures, require |E*| for a spectrum of temperature and frequency, i.e., require the definition of the master curve of the dynamic modulus. This is the case of two of the most internationally recognized design methods: CalME “California Mechanistic-Empirical Design” (3), and the “AASHTO Mechanistic-Empirical Pavement Design Guide”, known as the MEPDG (4).

The availability of dynamic modulus testing data corresponding to actual mixtures used in the pavement is not the most common situation. In many occasions, and particularly for secondary roads where the significance of the project does not justify additional testing expenses, only routine quality control test results are available. Even for highly important highway projects, actual mixtures are not available during the pavement design process. This led researchers to develop simplified models to predict dynamic modulus based on the percentage and the mechanical properties of the constituents of the asphalt mixtures. This approach has been introduced in the MEPDG for Levels 2 & 3 designs, when dynamic modulus data are not available but the composition of the actual mixture is known.

There are three widely accepted simplified predictive models. The first one is the Witczak dynamic modulus predictive equation (5), which is implemented in the MEPDG. It is a purely empirical model, which has been calibrated on the basis of 2750 data points from 205 different asphalt mixtures tested in the laboratories of the Asphalt Institute, the University of Maryland, and the Federal Highway Administration. Binder viscosity (η) is the rheological predictor variable in this model, which also includes aggregate gradation and mix volumetric parameters. The latter are voids in asphalt (VA) and percentage of voids in the mineral aggregate that are filled with asphalt (VFA). The second widely known and recognized simplified predictive model is the Hirsch model. This is a mechanical, though empirically-calibrated model, that is based on the law of mixtures (6). It combines the constituents of the asphalt mix (binder, voids, and aggregates) in a series and parallel phases arrangement. It was calibrated on the basis of 206 data points from 18 different asphalt mixtures tested at ALF-FHWA, MnROAD, and WesTrack full-scale testing facilities. Predictor variables in this model are VFA, VMA (voids in the mineral aggregate), and binder dynamic shear modulus (|G*|). The use of G* instead of η for binder rheology characterization was also adopted by Witczak and his coworkers in developing a second version of their predictive equation, which was calibrated on the basis of an extended database with 7400 data points from 346 different asphalt mixtures (7). This last model is usually referred to as Witczak 1-40D. It is worth to mention that the three volumetric factors in Witczak and Hirsch models (VA, VMA, and VFA) are routine mixture design parameters that are typically included in mixture design and specifications. Due to space limitations, actual formulas of the three models are not reproduced here.

1.1. Problem statementTOP

The empirical or semi-empirical nature of Witczak and Hirsch models recommends their validation for any database that differs from those used in the calibration process. A number of studies in the open literature present comparisons between these models and actual dynamic modulus data. Results from such studies do not always converge, with some of them reporting sufficiently accurate and robust estimates of |E*| for use in mechanistic-empirical pavement performance prediction and design, while others report the opposite. Examples of the former are provided by Schwartz (8), who evaluated Witczak η-based model and obtained an R² of 97% in the logarithmic scale, and Robbins (9), who reported Witczak η-based and Hirsch models predictions to generally follow the equality line when compared to measured data. Examples of unacceptable predictions for Witczak η-based model are provided by Birgisson et al. (10), who indicated the need to apply a correction factor that was between 0.668 and 2.740, depending on the method used for measuring binder viscosity, and Shen and Yu (11), who obtained modulus estimates that were between 3 and 7 times measured moduli for high temperatures and between 25% and 50% of measured moduli for low temperatures. Shen and Yu indicated that Witczak G*-based model improved predictions for high and low temperatures, while the opposite happened for intermediate temperatures. They also indicated that Hirsch model over predicts modulus at high temperatures. This over prediction has been also reported by Schwartz (8) for Witczak η-based model and by Al-Khateeb et al. (12) for both Witczak η-based and Hirsch models.

The above mentioned contradictory results, along with the fact that such a comparison effort has not yet been conducted for Spanish mixtures, requires caution when either of Witczak models or the Hirsch one is used for mixtures and testing protocols typically used in Spain. It should be mentioned that a number of the experimental studies presented above propose alternative simplified predictive equations as a result of Witczak and Hirsch models mispredictions (11, 12). Another problem is related to binder rheology characterization which, either in terms of η or G*, is not typically available in routine mixture design or specifications in Spain. The need exists to develop simplified predictive equations or a comprehensive database for such binder mechanical properties.

1.2. ObjectiveTOP

The objective of the experimental study presented herein is the evaluation of the applicability of Witczak and Hirsch models for asphalt concrete mixtures and testing protocols typically used in Spain. As a result of this evaluation, it is also considered the validation of an alternative simplified model that can be used to predict dynamic modulus of asphalt concrete from routine mixture design and specification parameters.

1.3. Research approachTOP

Dynamic modulus experimental data for the present research come from eight asphalt concrete mixtures (four different gradations and two asphalt contents) whose formulations are presented in Table 1. They represent typical mixtures used as base and binder courses in Spain. The number 20 in mixtures denomination denotes the nominal maximum aggregate size in mm. The letter “G” refers to a particular aggregate gradation that is only used as a base material, under a binder plus a wearing course. The letter “S” refers to another gradation that is typically used in Spain for binder or wearing courses. B60-70 is an unmodified asphalt binder whose present denomination, according to European normative EN 12591, is 50/70, which denotes the penetration interval (mm/10) at 25 °C. BM-3c (PMB 45/80-65, according EN 14023) is highly modified with SBS polymer, and is typically used in binder and wearing courses under heavy traffic. B13-22 (15/25, according EN 13924) is a hard-unmodified bitumen, which is used in the fabrication of high modulus asphalt concrete, “mezcla de alto módulo” (MAM) in Spanish.

Asphalt mixtures were produced in plant, and reheated in the laboratory in order to prepare gyratory specimens (ϕ100-mm; 180-mm height). Six specimens were prepared for each of the eight mixtures, and they were tested in stress-controlled cyclic compression in order to determine the complex modulus, E*. E* is a complex number whose modulus, |E*|, is the dynamic modulus and whose argument, φ, is referred to as phase angle, both being the so called viscoelastic properties under harmonic loading. Minimum testing temperature was −5 °C for all mixtures, while maximum temperature was adapted to binder rheology (50, 55, and 60 °C, respectively, for binders B60-70, BM-3c, and B13-22). A range of frequencies between 0.1 and 20 Hz were applied at temperature increments of no more than 15 °C. No rest period was introduced between two consecutive frequencies. Load level was adjusted for each frequency to produce ±35 με axial deformation, so that linear range was not exceeded. European standard EN 12697-26:2004 “Test Methods for Hot Mix Asphalt - Part 26: Stiffness”, annex D (cyclic tension-compression), was followed except for the type of loading, which was cyclic compression in the present study. The compressive cyclic loading used in this research is also specified by AASHTO T 342-11, whose methodology is very similar to that used in producing databases for Witczak and Hirsch models calibration. The only remarkable difference between the methodology used in the present experimental study and the AASHTO normative is the rest period between two consecutive frequencies, with the later specifying between 2 and 30 minutes. Coring of φ100-mm testing specimens from φ150-mm gyratory specimens is another prescription of AASHTO T 342 which is not specified by European normative and was not followed in this research.

As a first step, the measured dynamic moduli were compared to Witczak and Hirsch models predictions. This comparison involves 352 data points (44, on average, per mix), where each individual point is the mean of the six replicate specimens. Binder rheological properties to use in the models were determined according to AASHTO TP5-98 “Test Method for Determining the Rheological Properties of Asphalt Binder Using a Dynamic Shear Rheometer (DSR)”. The DSR device applies a torque to a thin film of asphalt binder placed between two parallel plates, while measuring the angular strain used to calculate the complex shear modulus, G*. A range of frequencies between 0.1 and 30 Hz were applied at temperatures between 4 and 82 °C, at 6 °C increments. It would have been desirable that the lower bound of the testing interval had been close to or below −5 °C, which was the lower bound of the corresponding temperature interval for asphalt mixture testing. Unfortunately, measuring binder stiffness at such low temperatures is beyond the loading capacity of the DSR used in this research.

2. CHARACTERIZATION OF MATERIALS IN THE LABORATORYTOP

2.1. Asphalt mixture complex modulusTOP

An example of dynamic moduli measured for one of the mixtures is presented in Figure 1. Values reflect temperature and frequency susceptibility of asphalt materials. Time-temperature superposition principle can be also deduced from the figure. According to this principle, which is applicable to thermorheologically simple materials like asphalt, a change in temperature is equivalent to a compression/extension of the time scale. This compression/extension means a shift in the logarithmic scale. Consequently, the sets of moduli measured for each temperature can be shifted in the log scale in order to produce a unique set, as reflected in Figure 1. This set defines a curve known as the “dynamic modulus master curve”, and the frequency resulting after the shifting process is known as the “reduced frequency”. The set was fitted by a sigmoidal function in this particular study. The same sigmoidal function, according to equation [1], is used by both CalME and the MEPDG.

where, |E*|  is  dynamic modulus

    f_red  is  reduced frequency: log (f_red) = log(f) +log (a_T)

    f   is  frequency (Hz)

    a_T  is  time-temperature shift factor; 20 °C is the reference temperature, for which a_T =1

    δ, α, β, & γ are model parameters

    note: log refers to decimal logarithm

Each of the parameters of the sigmoidal function has a distinct physical meaning. δ represents the minimum dynamic modulus of the mixture (10^δ), that is, the lower asymptote of the master curve. δ+α represents the maximum modulus (10^δ+α), which is the higher asymptote. β is a location parameter, which determines the horizontal position of the curvature sign change point, while γ represents the maximum slope of the master curve. Sigmoidal function is an analytical curve frequently used in mix/binder viscoelastic characterization within the frame of mechanistic-empirical pavement design. Its robust nature in curve fitting experimental data (both in frequency and time domain) with some noisy makes it also a good choice in pre-smoothing experimental data before fitting a more computational efficient analytical curve (i.e. Prony series) to perform stress-strain analysis of pavement layers using the Finite Element Method (13).

2.2. Binder rheologyTOP

Binder dynamic shear modulus data were also fitted by a master curve, after the application of the corresponding shift factors. An example is presented in Figure 2 for the B60-70 binder, where three sets of points can be observed. Each set corresponds to a different binder aging condition, which is mainly related to oxidation and volatile loss. The “original” condition refers to the binder before mixture fabrication; the “RTFO” condition refers to the binder residue after the artificially accelerated aging process of the Rolling Thin-Film Oven test (AASHTO T 240), conceived to reproduce aging during mixture fabrication; finally, the “PAV” condition refers to the residue after the Pressurized Aging Vessel test (AASHTO R 28), conceived to reproduce long-term aging (5–10 years) under service conditions. Binder mechanical properties used in mixture dynamic modulus predictive models must reflect actual binder aging condition in the mixture. In the present study, mixtures were tested after plant fabrication, so the binder aging condition is the one simulated by the RTFO test. CAM model (14), reflected in equations [2], was used to fit the |G*| master curve. The same model parameters were used for the three aging conditions, although a factor “k” was introduced that multiplies reduced frequency for RTFO and PAV data. This approach is based on aging effect being equivalent to a frequency shift in the logarithmic scale. The validity of this approach is based on the excellent fit that was achieved for original and PAV data while using a unique set of CAM model parameters, as observed in Figure 2. The same agreement was achieved for the rest of the binders.

where, |G*|  is  dynamic shear modulus

    δ   is  phase angle

    G_g  is  binder glassy modulus (1 GPa)

    f_red  is  reduced frequency

    R, w, & f_c are model parameters

Time-temperature shift factors obtained for G20 mixture and B60-70 binder are presented, as an example, in Figure 3. Agreement between both experimental sources was expected, and it happened for all temperatures except for the highest one. No definitive explanation was found for this disagreement, that was also observed for the rest of the mixtures except for those using B13-22, the hardest binder. This led authors to believe that the disagreement obtained for the highest temperature might be related to viscoplastic behavior, which is not captured by dynamic modulus. It is important to bear in mind that dynamic modulus can be considered a mechanical property of the asphalt mixture for medium and low temperatures, when linear viscoelasticity hypotheses are reasonably fulfilled. For high temperatures, non-linear viscoplastic behavior is expected, which means that the ratio of peak cyclic stress to peak cyclic strain will not be a constant anymore. Just as an example, Huang et al. (15) showed that confining pressure had a minimal effect on asphalt mixtures tested at 10 °C (linear behavior), while the opposite happened at 54 °C. At the last temperature, dynamic modulus values measured with confining pressures of 100 and 200 kPa were significantly greater than values measured without confining pressure, which indicates that dynamic modulus was not a material property.

Binder DSR minimum testing temperature was 5 °C, so low-temperature shift factors were extrapolated using the well-known Williams-Landel-Ferry equation reflected in Figure 3. The −5 °C factors obtained for the binders this way were in good agreement with values obtained for the asphalt mixtures. This result, together with the fact that all binders tend to reach a limiting modulus of 1 GPa for low temperatures, seem to overcome the DSR testing limitation.

Witczak first predictive equation is based on binder viscosity. This variable was determined by using the approach proposed by Bonaquist et al. (16), presented in equation [3], which is based on the classical Cox-Merz rule. It should be mentioned that binder behaves like a Newtonian fluid only at very high temperatures (17). For pavement service temperatures, result from equation [3] must be considered only an indication of binder stiffness rather than an intrinsic mechanical property.

where, η  is  binder viscosity in Pa·s

    |G*| is  dynamic shear modulus in Pa

    ω  is  angular frequency (rad/s)

    δ  is  phase angle

3. EVALUATION OF EXISTING PREDICTIVE EQUATIONSTOP

Both Witczak η-based and Hirsch models provided reasonable results when dynamic modulus predictions were compared to actual measured data. R² of Witczak η-based model was 96.2% in the logarithmic scale (88.9% in the arithmetic scale) and mean square error of the prediction (MSE) was 21.1% in the arithmetic scale. In the case of the Hirsch model, R² in the log scale was 95.7% (94.9% in the arithmetic scale) and MSE was 20.6%. These statistics are not very different from those reported by the authors of the two models in the calibration process (5, 6). Nonetheless, a clear pattern was observed in the prediction error of both Witczak η-based and Hirsch models. For high reduced frequencies (low temperatures), the former overestimated stiffness while the latter provided excellent predictions. For low reduced frequencies (high temperatures), both models underestimated moduli. This can be observed in the example presented in Figure 4, where measured and predicted moduli are plotted versus reduced frequency. The same pattern was observed for all mixtures, but not so distinctly for the high modulus asphalt concrete. This misprediction is related to differences between binder and mixture a_T shift factors for the highest temperature, that were explained above. It should be remarked that the high temperature pattern of the prediction error is not in line with a number of references that report both Witczak and Hirsch models to overpredict moduli under this condition (8, 11, 12).

Concerning Witczak 1-40D model, it clearly overestimated measured moduli. On average, predictions were 95.2% greater than measured values. R² was 76.3% in the log scale and MSE was 84.3% in the arithmetic scale. The same overprediction is reported by Robbins (9).

Another issue that can be deduced from Figure 4 is the low sensitivity of Witczak and Hirsch models versus binder content, since predictions are almost the same for both G20 and G20+ mixtures. The low sensitivity was observed for all the mixtures investigated in this study. This result was somewhat expected, since a number of studies exist where these models are shown to understate the influences of the mixture volumetric and gradation parameters (8, 11, 18). This fact, together with the pattern observed in error predictions, led authors to propose a new predictive equation where these two problems could be overcome.

4. PROPOSAL OF A NEW PREDICTIVE EQUATIONTOP

4.1. Equation calibrationTOP

The approach for the formulation of the new model is based on the relationships existing between parameters of the sigmoidal master curve (equation [1]) and mixture volumetric and gradation properties as well as binder rheological variables. There is a general agreement among researchers that δ and α parameters of the sigmoidal function are determined by volumetric and gradation properties of the mixtures rather than by binder rheology (8). In the upper limit of the reduced frequency (low temperatures), all binders become glassy materials with a modulus around 1 GPa, i.e., binder type does not make a difference. In the lower limit (high temperatures), binder stiffness becomes so low that aggregate skeleton determines asphalt mixture stiffness. Nonetheless, the efforts to relate δ and α parameters to mixture volumetric and gradation variables have not been very successful to date. The only variable for which consensus exists in the way it influences dynamic modulus is VMA, the voids in the mineral aggregate skeleton (6, 8, 12), which is included in Hirsch model but not in Witczak models. Unfortunately, the range of this variable was very narrow in this study, and such influence could not be directly determined. In order to overcome this limitation, results from this study were combined with those obtained by Shen and Yu (11) from seven mixtures typically used in Washington State. The sum of δ and α (log of maximum mixture modulus) is presented in Figure 5 for both experimental sources. The best fit equation presented in the figure can be used to extrapolate results from this study to mixtures with higher VMA. Extrapolation should not go much further than the actual VMA range tested in this research, which should be enough for most AC mixtures used in Spain. From Figure 5, Hirsch model understates the influence of VMA on dynamic modulus, which is in line with results presented above.

The sum δ+α varied from 4.43 to 4.55 among the different mixtures of this study, as reflected in Figure 5. In this case, the master curve of each mixture was determined independently from the others. The global R² of the master curve calibration was 99.93% in the logarithmic scale, which indicates an almost perfect fit, as can be seen in Figure 1. A second calibration was conducted where δ and α were forced to be equal for all the mixtures. The global R² of the new master curve calibration was still 99.82%, indicating that a single set of these two parameters can be used for all mixtures without losing accuracy. β parameter was found to depend on binder stiffness at reference temperature, whether it was formulated in terms of |G*| or η. This result was expected, as it has been pointed out by Schwartz (8) and reflected by Witczak models. But β was also found to depend on binder content, whether it was formulated in terms of VFA, VA, or effective binder volume (VMA-VA). This dependency is not acknowledged by Witczak models, and may be a reason for the low sensitivity of their predictions with respect to mixture volumetric properties. Finally, ϒ parameter was found to depend on β. Again, this dependency is not acknowledged by Witczak η-based model, which makes ϒ= −0.31335 in any case. The pattern described above was used to formulate a new predictive model that was calibrated from dynamic modulus data measured for the eight mixtures included in Table 1. The model is presented in equations [4], and Figure 6 indicates the model’s goodness of fit. R² of the model was 98.7% in the log scale (96% in arithmetic scale), and MSE was 12.6% in the arithmetic scale for the model based on G* and 12.5% for the model based on η.

δ+α=4.466

β= −1.177+9.316∙VA−0.5345∙log(G_ref) or +1.609+9.316∙VA−0.3352∙log(η_ref)

ϒ= −0.2094+0.09201∙β

log(a_T)= −0.1304∙(T−T_ref)

where, G_ref  is  binder dynamic shear modulus (MPa) at 20 °C (T_ref) and 10 rad/s

    η_ref  is  binder steady state viscosity (cP) at 20 °C (T_ref)

    VA  is  voids in asphalt (expressed as ratio)

    note: log refers to decimal logarithm

One limitation of the suggested model was related to using only three binder types in the calibration. Specifically, temperature susceptibility was relatively reduced in the three cases, with VTS (viscosity temperature susceptibility) parameter ranging from −2.85 to −2.75, which is a very narrow interval. VTS is the slope of the well-known relationship that exists between double logarithm of viscosity (cP) and logarithm of temperature (in Rankine). In order to overcome this limitation, an additional term was calibrated on the basis of a database generated by Onofre (19) at Universidade Federal do Ceará. Two asphalt binders with seven different types of modifiers were used by Onofre in his research. This way, the VTS varied from −4.25 to −2.90. As expected, temperature susceptibility influenced the ϒ parameter of the master curve, with higher temperature susceptibility corresponding to higher ϒ in absolute value. This influence could be modeled by adding the term “−0.1536·(|VTS|−2.824)” to the ϒ expression in equations [4]. This term can be used to extrapolate results from this study to mixtures with binders whose VTS parameter is below −2.85.

4.2. Practical implementation of the predictive equationTOP

Binder rheological data required by the new equation is very limited when compared to Witczak and Hirsch models, since binder stiffness must be only defined for a single temperature. But still, this information will not be available in daily practice of pavement design in Spain, since consistency prescriptions for binders are based on classical penetration and ring and ball tests rather than on DSR or viscosity testing.

A simplified approach is here presented in order to estimate η_ref and VTS from available data. First, it is assumed that the ring and ball softening point corresponds to a viscosity of 13000 Poises, which is a hypothesis based on research conducted by Shell Oil and later confirmed by Mirza and Witczak (20). Secondly, equation [5] is used to estimate log(η_ref) from binder penetration. This empirical equation was calibrated (R²=97.5%; MSE=1.17%) on the basis of ten binders typically used in Spain, that were tested within the frame of a research project focused on binder rheology (21). Once viscosity is known for two different temperatures, the VTS parameter can be determined.

where, η_ref is  binder steady state viscosity (cP) at 20 °C

    P  is  binder penetration (mm/10) at 25 °C

    note: log refers to decimal logarithm

4.3. Validation of the predictive equationTOP

Predictions of the new equation were compared to dynamic modulus data measured for five mixtures coming from the CEDEX accelerated pavement testing facility. Mixtures had been tested between −15 and 40 °C. The simplified approach described above was used to estimate η_ref and VTS from binders penetration and ring and ball softening point. The RTFO test residue data were used, since asphalt mixtures were tested after fabrication (no additional aging had taken place). Model’s goodness of fit is appreciated in Figure 7. R² of the prediction was 97.9% in the log scale (94.1% in the arithmetic scale) and MSE was 14.4% in the arithmetic scale. No bias is appreciated in the error. These statistics are similar to those obtained in the calibration, which validates the suggested predictive equation.

5. CONCLUSIONSTOP

Dynamic modulus has been measured for eight asphalt concrete mixtures that were tested between −5 and 60 °C. They represent typical mixtures used as base and binder courses in Spain. These experimental data have been used to evaluate three widely accepted dynamic modulus predictive models: Witczak η-based equation, which is formulated in terms of binder viscosity, and Witczak 1-40D and Hirsch models, which are based on binder dynamic shear modulus. Witczak 1-40D equation considerably overestimated (95% on average) measured dynamic modulus. On the contrary, both Witczak η-based and Hirsch models provided reasonable results. For both models, R² was 96% in the log scale and the mean square error of the prediction was 21% in the arithmetic scale, which are not very different from the statistics reported by the authors in the calibration process. Nonetheless, a clear pattern was observed in the prediction error, with both models underestimating modulus for low reduced frequencies (high temperatures) while Witczak η-based model overestimating modulus for high reduced frequencies (low temperatures). Besides, the three models presented very low sensitivity versus binder content, which was against experimental results.

A new predictive model was formulated in order to overcome the limitations mentioned above. The suggested model is based on the relationships existing between the four parameters of the sigmoidal master curve (δ, α, β, and ϒ) and mixture volumetric and gradation properties as well as binder rheological variables. The model was calibrated on the basis of 352 data points from the eight asphalt concrete mixtures. R² of the calibration process was 98.7% in the log scale (96% in arithmetic scale), and the mean square error of the prediction was 12.5% in the arithmetic scale.

Binder rheological data required by the new model is very limited when compared to Witczak and Hirsch models, since binder stiffness must be only defined for a single temperature, 20 °C, either in terms of the steady state viscosity or in terms the dynamic shear modulus. Binder viscosity temperature susceptibility parameter (VTS) must be also defined. The other predictor variables are voids in asphalt, and voids in the mineral aggregate, which are routine mixture design parameters. Besides, a simplified procedure is proposed for estimating such binder data from penetration and ring and ball test results.

The new model was validated, together with the simplified procedure for binder data estimation, on the basis of the dynamic modulus values measured for five asphalt concrete mixtures that had been tested, within the frame of other research projects, between −15 and 40 °C. R² of the prediction was 97.9% in the log scale and the mean square error was 14.4% in the arithmetic scale. No bias was found in the error.

ACKNOWLEDGEMENTSTOP

Authors would like to express their gratitude for the important laboratory effort behind this study. Workers who made it possible are sincerely acknowledged and recognized.

REFERENCESTOP