Development and Analysis of High-Modulus Asphalt Concrete Predictive Model

The main purpose of this paper is to present the development of a new predictive model intended for the calculation of stiffness modulus |E*| determined by a four-point bending beam test (4PBB or 4PB-PR). The model developed, called model A, was based on the Witczak model, which was developed for the dynamic-modulus (DM) method. Most of the asphalt mixtures used to develop the model were high-modulus asphalt concrete (HMAC). The most commonly used methods for determining the stiffness modulus |E*| of asphalt mixtures were also discussed. The paper presents the results of the study for 10 asphalt mixtures but 8 of them were used to develop the predictive model. In addition, the results of complex shear modulus G* tests on neat and modified bituminous binders carried out in a dynamic shear rheometer (DSR), necessary for the development of a predictive model, are presented. The tests carried out in the dynamic shear rheometer had significant measurement uncertainties. The results of the volumetric parameters of the asphalt mixtures are also reported. The developed model A has maximum absolute errors e = 1930 MPa (p = 95%) and maximum relative errors re = 50% (p = 95%). The distribution of the absolute errors of the model, after discarding outliers, has a normal distribution as in the development of other models of this type, which was confirmed by appropriate statistical tests. On the basis of the tests and calculations carried out, it was concluded that, in order to increase the precision of the predictive models, it is advisable to reduce the measurement uncertainty of the bitumen complex shear modulus G*. For the developed model A, the limiting values of the stiffness modulus |E*| are also shown, within which the determined stiffness modulus should fall.


Introduction
The stiffness modulus |E*| of an asphalt mixture is an important parameter characterizing its properties. Different values of stiffness modulus |E*| appear at particular stages of asphalt pavement design and realization. The stiffness modulus |E*| is used in the designing process of pavements through mechanistic methods. It is also determined in the case of asphalt mixture composition designing as well as the factory production-control process [1].
The stiffness modulus |E*| is a physical property describing the relationship between stress and strain in the case of linear viscoelastic material as it is loaded and unloaded. When the stress and stiffness modulus are known, the strain and displacement can be calculated. The stiffness modulus |E*| extensively describes the asphalt mix, due to its dependence on many factors. These factors are the composition of an asphalt mixture, the content of air voids, shape, lay down, and composition of aggregate grains in an asphalt mixture. The values of the stiffness modulus |E*| change substantially along with the temperature and loading frequency [2,3] which makes the cause of differentiation between the stiffness modulus |E*| and modulus of elasticity. The value of stiffness modulus |E*| depends on the method of determining, testing equipment, conditions during the tests, and 1.
Determination of stiffness modulus and phase angle of an asphalt mixture at various temperatures and various loading frequencies for each test temperature by the fourpoint bending beam test conducted according to [6]; 2.
Determination of shear modulus and phase angle of a binder at various temperatures and various loading frequencies for each test temperature by a dynamic shear rheometer test, carried out according to [15]; 3.
Determination of the resistance to hardening of binders under the influence of heat and air by an RTFOT test performed according to [16]; 4.
Determination of soluble binder content in an asphalt mixture conducted according to [17]; 5.
Determination of particle size distribution carried out according to [18]; 6.
Determination of the maximum density of an asphalt mixture according to [19]; 7.
Determination of bulk density of an asphalt mixture specimen according to [20]; 8.
Determination of air-voids content of an asphalt-mixture specimen according to [21]; 9.
Determination of dimensions of an asphalt-mixture specimen according to [22].
During the tests, every effort was made to meet the requirements of the abovementioned standards.

Tested Asphalt Mixtures and Binders
Ten asphalt mixtures were tested as part of the research project. Three HMA (hotmix asphalt) types were developed and tested specifically for the project, the other seven ones were tested as part of the standard activities of the Road Laboratory of the General Directorate of National Roads and Motorways in Poland. The testing program of the standard-tested HMA for the research project was extended. Most of the relevant data for the tested HMA are presented in Table 1. HMA mixtures were made with both unmodified and polymer-modified binders. Many research works concerning polymer-modified binders conducted previously by the authors were described in [23][24][25]. The gradation curves of aggregate mixtures used in each HMA are shown in Figure 1. The mixture with the lowest voids content is plotted in green (HMA 10) while the mixture with the highest voids content is plotted in red (HMA 6). The limits for the gradation curves are given in accordance with the Polish technical requirements for high-modulus asphalt concrete [26].
Among the ten HMA types tested, taking into account the Polish requirements [26], seven mixtures can be considered HMAC; these are mixtures in which mainly unmodified hard bitumen 20/30, or in one case with polymer modified bitumen PMB 25/55-60, were used. HMA 8 and HMA 10 mixtures contain a bituminous binder suitable for HMAC and the gradation curve also met the requirements specified for HMAC [26] but the values of stiffness modulus determined at 10 • C and 10 Hz were found to be less than E 4PB = 11,000 MPa, the minimum stiffness modulus value specified for HMAC for the base course, and, therefore, HMA 8 and HMA 10 mixtures, were recognized as asphalt concrete (AC). HMA 9, which is an SMA Jena mix (stone mastic asphalt used as a single-layer pavement), was decided to be included in the research program as a comparative mix. The six HMA mixtures (HMA 1-HMA 6) contain different 20/30 penetration-grade bitumen from three different manufacturers. HMA 7 and HMA 9 mixtures contain polymer-modified bitumen and HMA 8 and HMA 10 mixtures used PMB 25/55-80 highly modified bitumen (so-called HIMA), both from the same manufacturer. In all the tested HMA mixtures, limestone filler was used. An adhesion agent was added to all the tested HMA in an amount ranging from 0.3% to 0.5% by weight of the bitumen. Apart from the adhesive agent, no other modifiers were added to the HMA [27].  Among the ten HMA types tested, taking into account the Polish requirements [26], seven mixtures can be considered HMAC; these are mixtures in which mainly unmodified hard bitumen 20/30, or in one case with polymer modified bitumen PMB 25/55-60, were used. HMA 8 and HMA 10 mixtures contain a bituminous binder suitable for HMAC and the gradation curve also met the requirements specified for HMAC [26] but the values of stiffness modulus determined at 10 °C and 10 Hz were found to be less than E4PB = 11,000 MPa, the minimum stiffness modulus value specified for HMAC for the base course, and, therefore, HMA 8 and HMA 10 mixtures, were recognized as asphalt concrete (AC). HMA 9, which is an SMA Jena mix (stone mastic asphalt used as a single-layer pavement), was decided to be included in the research program as a comparative mix. The six HMA mixtures (HMA 1-HMA 6) contain different 20/30 penetration-grade bitumen from three different manufacturers. HMA 7 and HMA 9 mixtures contain polymer-modified bitumen and HMA 8 and HMA 10 mixtures used PMB 25/55-80 highly modified bitumen (so-called HIMA), both from the same manufacturer. In all the tested HMA mixtures, limestone filler was used. An adhesion agent was added to all the tested HMA in an amount ranging from 0.3% to 0.5% by weight of the bitumen. Apart from the adhesive agent, no other modifiers were added to the HMA [27].
Finally, the results obtained for eight HMA mixtures (HMA 1-HMA 8) were used to develop the predictive model, of which seven were classified as HMAC and one (HMA 8) as asphalt concrete. Mixtures HMA 9 and HMA 10 were not included in the development of the predictive model, as their properties differed too much from the requirements specified for HMAC.
The asphalt mixtures used in the study were prepared in a laboratory mixer ( Figure  2), with the exception of HMA 6, which was taken from the site and transported to the Finally, the results obtained for eight HMA mixtures (HMA 1-HMA 8) were used to develop the predictive model, of which seven were classified as HMAC and one (HMA 8) as asphalt concrete. Mixtures HMA 9 and HMA 10 were not included in the development of the predictive model, as their properties differed too much from the requirements specified for HMAC.
The asphalt mixtures used in the study were prepared in a laboratory mixer ( Figure 2), with the exception of HMA 6, which was taken from the site and transported to the laboratory. The mixing and paving temperatures were adopted depending on the type of mix and the type of bitumen used according to [26]. Two slabs of HMA from each mix were compacted in a roller device (Figure 3a) according to [28]. Four prismatic specimens ( Figure 3b) with dimensions 50 mm × 60 mm × 380 mm (H × B × L) were cut from the compacted slabs, on which the stiffness modulus and phase angle could be determined using the four-point bending test (4PB-PR).
Each bitumen applied in the tested HMA was subjected to short-term RTFOT ageing in accordance with [16] prior to testing in the DSR.  Each bitumen applied in the tested HMA was subjected to short-term RTFOT ageing in accordance with [16] prior to testing in the DSR.

Shear Modulus and Phase Angle of Bituminous Binders
The shear modulus and phase angle of the investigated bituminous binders were determined using a dynamic shear rheometer DSR, shown in Figure 4. Each type of bitumen used for the analyzed HMA was tested in the DSR. Tests were carried out at different temperatures and loading frequencies for each temperature so that the temperature and loading frequency on the DSR device were consistent with the temperature and loading frequency used in the determination of the stiffness modulus and phase angle of the HMA mixtures. The selection of test temperatures also took into account the limiting operating (equivalent) temperatures according to the Superpave classification. The temperatures and loading frequencies of the individual binders are shown in Table 2. In the first series, the tests were also carried out at high temperatures using a 25 mm diameter spindle. However, after analyzing the first series of tests, it became apparent that some of the test results may not meet the measurement uncertainty requirements specified in [15]. Therefore, it was decided to amend the test procedure and carry out a second series of tests. The second series of tests was carried out using only the 8 mm diameter spindle in the temperature range of −20°C to 40 °C.  Each bitumen applied in the tested HMA was subjected to short-term RTFOT ageing in accordance with [16] prior to testing in the DSR.

Shear Modulus and Phase Angle of Bituminous Binders
The shear modulus and phase angle of the investigated bituminous binders were determined using a dynamic shear rheometer DSR, shown in Figure 4. Each type of bitumen used for the analyzed HMA was tested in the DSR. Tests were carried out at different temperatures and loading frequencies for each temperature so that the temperature and loading frequency on the DSR device were consistent with the temperature and loading frequency used in the determination of the stiffness modulus and phase angle of the HMA mixtures. The selection of test temperatures also took into account the limiting operating (equivalent) temperatures according to the Superpave classification. The temperatures and loading frequencies of the individual binders are shown in Table 2. In the first series, the tests were also carried out at high temperatures using a 25 mm diameter spindle. However, after analyzing the first series of tests, it became apparent that some of the test results may not meet the measurement uncertainty requirements specified in [15]. Therefore, it was decided to amend the test procedure and carry out a second series of tests. The second series of tests was carried out using only the 8 mm diameter spindle in the temperature range of −20°C to 40 °C.

Shear Modulus and Phase Angle of Bituminous Binders
The shear modulus and phase angle of the investigated bituminous binders were determined using a dynamic shear rheometer DSR, shown in Figure 4. Each type of bitumen used for the analyzed HMA was tested in the DSR. Tests were carried out at different temperatures and loading frequencies for each temperature so that the temperature and loading frequency on the DSR device were consistent with the temperature and loading frequency used in the determination of the stiffness modulus and phase angle of the HMA mixtures. The selection of test temperatures also took into account the limiting operating (equivalent) temperatures according to the Superpave classification. The temperatures and loading frequencies of the individual binders are shown in Table 2. In the first series, the tests were also carried out at high temperatures using a 25 mm diameter spindle. However, after analyzing the first series of tests, it became apparent that some of the test results may not meet the measurement uncertainty requirements specified in [15]. Therefore, it was decided to amend the test procedure and carry out a second series of tests. The second series of tests was carried out using only the 8 mm diameter spindle in the temperature range of −20 • C to 40 • C.
swing angle value of θ = 0.001 rad in the frequency range of 0.1 Hz to 50 Hz. Tests were conducted from lowest to highest frequency. The temperature gradient was not greater than 1 °C/min. For each type of binder, the tests were carried out on four samples with a diameter of 8 mm and on four samples with a diameter of 25 mm. As a result of the test, the measurement system was saved to a .csv file on the control computer. Figures 5-9 show chosen test results concerning shear modulus |G*| and phase angle δ of asphalt binders at loading frequency 1.59 Hz (10 rad/s).   Tests at each temperature and loading frequency were carried out on four bitumen samples. After a detailed review of the standards [15,29,30] for testing in the DSR apparatus, it was decided that the determinations of shear modulus and phase angle would be carried out from the highest to the lowest temperature. It was considered that the extreme test conditions, i.e., highest and lowest temperature and lowest and highest frequency, would be the most dangerous for the specimen which, after testing, was reflected in a higher measurement uncertainty for such conditions. When testing at high temperature, it was observed whether the specimen (at constant strain amplitude and maximum frequency) was still within the range of linear strain. Liquefaction and spillage of the sample on the lower measuring plate of the instrument resulted in the rejection of the test results. When conducting low-temperature tests, it was observed whether the sample still had proper adhesion to the lower measuring plate after the test. If it was found that adhesion had been lost, the results were discarded. Tests were carried out at a constant maximum swing angle value of θ = 0.001 rad in the frequency range of 0.1 Hz to 50 Hz. Tests were conducted from lowest to highest frequency. The temperature gradient was not greater than 1 • C/min. For each type of binder, the tests were carried out on four samples with a diameter of 8 mm and on four samples with a diameter of 25 mm. As a result of the test, the measurement system was saved to a .csv file on the control computer. Figures 5-9 show chosen test results concerning shear modulus |G*| and phase angle δ of asphalt binders at loading frequency 1.59 Hz (10 rad/s).
During the preparation for the implementation of the study, the authors obtained information about the low reproducibility of test results obtained from DSR. Such news was also confirmed by the staff of the main laboratory of one of the companies. The data on the precision of the MSCR test available in the standard [30], i.e., a reproducibility of up to 43%, also did not inspire optimism. After carrying out the first series of tests, it turned out that the obtained results of the shear modulus |G*| were characterized by too high of a measurement uncertainty, significantly exceeding the 10% specified in [15] and did not meet the requirement for the arithmetic mean at the overlap temperature. It was decided to introduce changes to the DSR test procedure for the binders. Particular attention was paid to maintaining the correct specimen geometry and ensuring good adhesion of the binder to the rheometer test system by heating both test-system plates.
These problems are mentioned in the European standard [15] but no specific solutions to these problems are given. After modifying the measurement procedure, a second batch of binder tests was carried out, with significantly less measurement uncertainty. The results obtained were assessed according to the reproducibility criteria contained in [15]. The criteria were considered to be fulfilled if the mean relative uncertainty RMU of the results of the determination of the shear modulus |G*| and the phase angle δ was less than 10% (for the shear modulus |G*|) and 5% (for the phase angle δ), respectively. The largest values of relative measurement uncertainty are obtained for extreme temperature and frequency conditions.  During the preparation for the implementation of the study, the authors obtained information about the low reproducibility of test results obtained from DSR. Such news was also confirmed by the staff of the main laboratory of one of the companies. The data on the precision of the MSCR test available in the standard [30], i.e., a reproducibility of up to 43%, also did not inspire optimism. After carrying out the first series of tests, it turned out that the obtained results of the shear modulus |G*| were characterized by too high of a measurement uncertainty, significantly exceeding the 10% specified in [15] and did not meet the requirement for the arithmetic mean at the overlap temperature. It was decided to introduce changes to the DSR test procedure for the binders. Particular attention was paid to maintaining the correct specimen geometry and ensuring good adhesion of the binder to the rheometer test system by heating both test-system plates. These problems are mentioned in the European standard [15] but no specific solutions to these  During the preparation for the implementation of the study, the authors obtained information about the low reproducibility of test results obtained from DSR. Such news was also confirmed by the staff of the main laboratory of one of the companies. The data on the precision of the MSCR test available in the standard [30], i.e., a reproducibility of up to 43%, also did not inspire optimism. After carrying out the first series of tests, it turned out that the obtained results of the shear modulus |G*| were characterized by too high of a measurement uncertainty, significantly exceeding the 10% specified in [15] and did not meet the requirement for the arithmetic mean at the overlap temperature. It was decided to introduce changes to the DSR test procedure for the binders. Particular attention was paid to maintaining the correct specimen geometry and ensuring good adhesion of the binder to the rheometer test system by heating both test-system plates. These problems are mentioned in the European standard [15] but no specific solutions to these criteria were considered to be fulfilled if the mean relative uncertainty RMU of the results of the determination of the shear modulus |G*| and the phase angle δ was less than 10% (for the shear modulus |G*|) and 5% (for the phase angle δ), respectively. The largest values of relative measurement uncertainty are obtained for extreme temperature and frequency conditions.   criteria were considered to be fulfilled if the mean relative uncertainty RMU of the results of the determination of the shear modulus |G*| and the phase angle δ was less than 10% (for the shear modulus |G*|) and 5% (for the phase angle δ), respectively. The largest values of relative measurement uncertainty are obtained for extreme temperature and frequency conditions.  It is advisable to develop measurement-precision data for the entire spectrum of temperatures and frequencies (0.1 to 25 Hz), with a particular focus on frequencies of 1.59 Hz and 10 Hz.
The results were compared for the five selected temperatures used in the tests: −20 • C, 0 • C, 10 • C, 40 • C, and 60 • C. At low and medium temperatures, the results of the determined shear modulus can be divided in terms of obtained values into three groups:  It is advisable to develop measurement-precision data for the entire spectrum of temperatures and frequencies (0.1 to 25 Hz), with a particular focus on frequencies of 1.59 Hz and 10 Hz.
The results were compared for the five selected temperatures used in the tests: −20 °C, 0 °C, 10 °C, 40  At 0 °C, the difference between the shear modulus determined for the B1 binder (111.8 MPa) and the average determined for the rest of the hard binders (69.7 MPa) is 42 MPa or 60% of the lower value. At high temperatures, the difference between the B1 binder and the other 20/30 unmodified binders disappears. However, it is still possible to observe clear differences in the values obtained for neat and polymer-modified binders. The shear modulus values decrease with increasing temperature and increase with increasing frequency, as can be observed in other studies. The obtained values confirm the fact that modified binders have a lower temperature susceptibility.
In the case of the phase angle of the tested binders, also at low and medium temperatures, a clear difference can be seen between the 20/30 unmodified binders with smaller phase angle values and the modified binders. At high temperature, on the other hand, it is the 20/30 unmodified binders that are characterized by larger phase-angle values. At 60 °C, the average value of the phase angle for hard binders (65.9°) is 8° higher than the average for polymer-modified binders (57.9°). The values of the phase angle of binders increase with increasing temperature and decrease with increasing frequency. The smallest value of the phase angle (16.2°) was obtained for B1 bitumen at −20 °C. The largest value of phase angle (68.4°) was also obtained for B1 bitumen at 60 °C.

Stiffness Modulus and Phase Angle of Asphalt Mixtures
Determination of the stiffness modulus and phase angle of the HMA mixtures was carried out by using the 4PBB bending beam apparatus ( Figure 10). The tests were performed at different temperatures and frequencies for each temperature. The variation in At 0 • C, the difference between the shear modulus determined for the B1 binder (111.8 MPa) and the average determined for the rest of the hard binders (69.7 MPa) is 42 MPa or 60% of the lower value. At high temperatures, the difference between the B1 binder and the other 20/30 unmodified binders disappears. However, it is still possible to observe clear differences in the values obtained for neat and polymer-modified binders. The shear modulus values decrease with increasing temperature and increase with increasing frequency, as can be observed in other studies. The obtained values confirm the fact that modified binders have a lower temperature susceptibility.
In the case of the phase angle of the tested binders, also at low and medium temperatures, a clear difference can be seen between the 20/30 unmodified binders with smaller phase angle values and the modified binders. At high temperature, on the other hand, it is the 20/30 unmodified binders that are characterized by larger phase-angle values. At 60 • C, the average value of the phase angle for hard binders (65.9 • ) is 8 • higher than the average for polymer-modified binders (57.9 • ). The values of the phase angle of binders increase with increasing temperature and decrease with increasing frequency. The smallest value of the phase angle (16.2 • ) was obtained for B1 bitumen at −20 • C. The largest value of phase angle (68.4 • ) was also obtained for B1 bitumen at 60 • C.

Stiffness Modulus and Phase Angle of Asphalt Mixtures
Determination of the stiffness modulus and phase angle of the HMA mixtures was carried out by using the 4PBB bending beam apparatus ( Figure 10). The tests were performed at different temperatures and frequencies for each temperature. The variation in test temperature is intended to represent the different operating temperatures of the pavement. The different test frequencies also mimic the frequencies found in the pavement, which are primarily dependent on traffic speed.
The tests were carried out on 40 HMA beam samples made from 10 HMA mixtures (four samples per one HMA). All beam samples had the same dimensions, i.e., 50 mm × 60 mm × 380 mm (H × B × L). The specimens were conditioned from 4 to 6 h before the start of the test at each temperature. The execution of tests at different frequencies could be set in the test apparatus control program. The tests were carried out in controlled-strain mode, which in all cases was ε = 50 µm/m, thus avoiding fatigue processes. For each determination, 150 loading cycles were performed. Fifteen cycles (from 93 to 107 cycles) were used to calculate the values of the stiffness modulus E 4PB and the phase angle Φ 4PB of the asphalt mixture. test temperature is intended to represent the different operating temperatures of the pavement. The different test frequencies also mimic the frequencies found in the pavement, which are primarily dependent on traffic speed. The tests were carried out on 40 HMA beam samples made from 10 HMA mixtures (four samples per one HMA). All beam samples had the same dimensions, i.e., 50 mm × 60 mm × 380 mm (H × B × L). The specimens were conditioned from 4 to 6 h before the start of the test at each temperature. The execution of tests at different frequencies could be set in the test apparatus control program. The tests were carried out in controlled-strain mode, which in all cases was ε = 50 μm/m, thus avoiding fatigue processes. For each determination, 150 loading cycles were performed. Fifteen cycles (from 93 to 107 cycles) were used to calculate the values of the stiffness modulus E4PB and the phase angle Φ4PB of the asphalt mixture.
Temperature and frequency were not the same for each mixture due to technical reasons. Table 3 summarizes the temperature and loading frequency used in the tests for each HMA. The lowest temperature at which the mixtures were tested was 0 °C, as the climatic chamber of the testing machine did not allow a lower temperature to be set reliably. The highest temperature at which the mixtures were tested was 40 °C. According to [6], this is the highest temperature at which specimens should be tested due to the possibility of nonlinear deformation and the creep phenomenon of HMA. The tests were carried out from the lowest to the highest temperature and from the lowest to the highest frequency. For HMA 1, only the standard tests required by [26] were made. Due to the higher probability of nonlinear deformations, HMA 9 samples were not tested at 40 °C.  Temperature and frequency were not the same for each mixture due to technical reasons. Table 3 summarizes the temperature and loading frequency used in the tests for each HMA. The lowest temperature at which the mixtures were tested was 0 • C, as the climatic chamber of the testing machine did not allow a lower temperature to be set reliably. The highest temperature at which the mixtures were tested was 40 • C. According to [6], this is the highest temperature at which specimens should be tested due to the possibility of nonlinear deformation and the creep phenomenon of HMA. The tests were carried out from the lowest to the highest temperature and from the lowest to the highest frequency. For HMA 1, only the standard tests required by [26] were made. Due to the higher probability of nonlinear deformations, HMA 9 samples were not tested at 40 • C. The tests discussed here yielded 855 determinations of the complex stiffness modulus E*, after excluding repetitive frequencies and files containing errors and the results obtained for HMA 9 and HMA 10 to develop an analytical-empirical model, 471 values of the complex stiffness modulus E* were obtained from the eight tested HMA. Figures 11-13 show chosen test results of stiffness modulus |E*| of hot-mix asphalt for frequency 10 Hz at different temperatures. Figure 12 shows the limiting values of stiffness modulus |E*| specified in the Polish technical requirements [26].
The tests discussed here yielded 855 determinations of the complex stiffness modulus E*, after excluding repetitive frequencies and files containing errors and the results obtained for HMA 9 and HMA 10 to develop an analytical-empirical model, 471 values of the complex stiffness modulus E* were obtained from the eight tested HMA. Figures  11-13 show chosen test results of stiffness modulus |E*| of hot-mix asphalt for frequency 10 Hz at different temperatures. Figure 12 shows the limiting values of stiffness modulus |E*| specified in the Polish technical requirements [26].   The tests discussed here yielded 855 determinations of the complex stiffness modulus E*, after excluding repetitive frequencies and files containing errors and the results obtained for HMA 9 and HMA 10 to develop an analytical-empirical model, 471 values of the complex stiffness modulus E* were obtained from the eight tested HMA. Figures  11-13 show chosen test results of stiffness modulus |E*| of hot-mix asphalt for frequency 10 Hz at different temperatures. Figure 12 shows the limiting values of stiffness modulus |E*| specified in the Polish technical requirements [26].   The measurement results (expressed by mean value with uncertainty interval) of the stiffness modulus and phase angle were determined on a sample of 15 values. The largest measurement uncertainty values (approximately 10 MPa) were found in the case of the tests carried out at the lowest applied frequency of 0.1 Hz. All mean relative measurement uncertainty values were less than 1%.
When the precision of the determination of the stiffness modulus and the phase angle using the four-point bending beam method (4PBB) is considered, it should be borne in mind that a significant source of uncertainty is the way the specimen is mounted in the test device. This was also mentioned in [31]. The value of the measurement uncertainty of these quantities is also influenced by the fact that the algorithm for approximation of the measured data is not specified in [6], and the possible different possibilities of calculating the result of the determination (using the approximation, or directly from the measured values) that result from this fact. The measurement results (expressed by mean value with uncertainty interval) of the stiffness modulus and phase angle were determined on a sample of 15 values. The largest measurement uncertainty values (approximately 10 MPa) were found in the case of the tests carried out at the lowest applied frequency of 0.1 Hz. All mean relative measurement uncertainty values were less than 1%.
When the precision of the determination of the stiffness modulus and the phase angle using the four-point bending beam method (4PBB) is considered, it should be borne in mind that a significant source of uncertainty is the way the specimen is mounted in the test device. This was also mentioned in [31]. The value of the measurement uncertainty of these quantities is also influenced by the fact that the algorithm for approximation of the measured data is not specified in [6], and the possible different possibilities of calculating the result of the determination (using the approximation, or directly from the measured values) that result from this fact.
At low temperature (0 °C), the higher value of the stiffness modulus E4PB = 21,080 MPa (mean) is achieved by HMA 7 and HMA 8 with a modified binder with 5.2% bitumen content by weight, while HMA 6 with an unmodified binder with the same bitumen content by weight has a significantly lower (by 3880 MPa) value of the stiffness modulus E4PB = 17,200 ± 548 MPa, which is, according to the authors' opinion, an unusual behaviour. At medium temperature (10 °C), the highest value of stiffness modulus E4PB = 18,530 ± 1137 MPa is achieved by HMA 1, which corresponds well with the value of shear modulus of B1 bitumen. According to statistical tests performed at 10 °C, the stiffness modulus E4PB of the HMA 1 mix has a higher value than that of HMA 3 (17,606 ± 907 MPa), HMA 4 (16,459 ± 1219 MPa), and HMA 5 (17,698 ± 613 MPa) mixes. In contrast, the E4PB stiffness modulus values of the HMA 3 and HMA 5 mixtures are equal to each other. However, the relative difference (5%) between the E4PB stiffness modulus values of HMA 1 and HMA 3 and HMA 5 is much smaller than the relative difference (75%) between the shear modulus values of the bitumen used as a binder, as B1 bitumen has a significantly higher shear modulus value |G*|. This can be explained by the higher volumetric content of bitumen in HMA 1 and the differences in voids content Va. At 40 °C, HMA 6 with 20/30, bitumen stands out with a significantly higher value of stiffness modulus E4PB = 2925 ± 429 MPa than mixtures with polymer-modified bitumens (average E4PB = 1239 MPa). At low temperature (0 • C), the higher value of the stiffness modulus E 4PB = 21,080 MPa (mean) is achieved by HMA 7 and HMA 8 with a modified binder with 5.2% bitumen content by weight, while HMA 6 with an unmodified binder with the same bitumen content by weight has a significantly lower (by 3880 MPa) value of the stiffness modulus E 4PB = 17,200 ± 548 MPa, which is, according to the authors' opinion, an unusual behaviour. At medium temperature (10 • C), the highest value of stiffness modulus E 4PB = 18,530 ± 1137 MPa is achieved by HMA 1, which corresponds well with the value of shear modulus of B1 bitumen. According to statistical tests performed at 10 • C, the stiffness modulus E 4PB of the HMA 1 mix has a higher value than that of HMA 3 (17,606 ± 907 MPa), HMA 4 (16,459 ± 1219 MPa), and HMA 5 (17,698 ± 613 MPa) mixes. In contrast, the E 4PB stiffness modulus values of the HMA 3 and HMA 5 mixtures are equal to each other. However, the relative difference (5%) between the E 4PB stiffness modulus values of HMA 1 and HMA 3 and HMA 5 is much smaller than the relative difference (75%) between the shear modulus values of the bitumen used as a binder, as B1 bitumen has a significantly higher shear modulus value |G*|. This can be explained by the higher volumetric content of bitumen in HMA 1 and the differences in voids content V a . At 40 • C, HMA 6 with 20/30, bitumen stands out with a significantly higher value of stiffness modulus E 4PB = 2925 ± 429 MPa than mixtures with polymer-modified bitumens (average E 4PB = 1239 MPa).

Effective Bitumen Content, Air Voids, and Granulation of HMA
One of the variables used in the equation of the Witczak-El-Badawy model is the effective bitumen content in the HMA V beff . In order to calculate the value of this variable, it is necessary to determine the density of bitumen G b used in the HMA. The authors adopted the values of bitumen density G b from data made available by bitumen manufacturers. The results of the effective bitumen content V beff obtained from tests carried out on the basis of [17] are shown in Figure 14. method B (in water) was used to calculate the air-voids content. The bulk density of specimens can also be determined using method D (from sample dimensions). Howev on the basis of the tests carried out, it can be assumed that due to the small difference values, the determination method has no influence on the bulk density value. The resu of the air-voids content obtained from the tests [20] using method B (in water) are p sented in Figure 15. This figure shows also the limiting values of air voids Va specified the Polish technical requirements [23].  From the sensitivity analysis performed according to [11], it appears that the value of the stiffness modulus depends on the air-voids content of the HMA. In view of the mixes with the lowest homogeneity, it was decided that the stiffness modulus values marked on the individual beams with their corresponding air-voids contents would be used to optimize the predictive-model coefficients. The bulk density determined by method B (in water) was used to calculate the air-voids content. The bulk density of the specimens can also be determined using method D (from sample dimensions). However, on the basis of the tests carried out, it can be assumed that due to the small difference in values, the determination method has no influence on the bulk density value. The results of the air-voids content obtained from the tests [20] using method B (in water) are presented in Figure 15. This figure shows also the limiting values of air voids Va specified in the Polish technical requirements [23].
The granulation data was calculated based on the gradation-curve values determined for the HMA used for the beam specimens ( Table 4). Explanations of the grain-size data are given next to Equation (2).  The granulation data was calculated based on the gradation-curve values determined for the HMA used for the beam specimens ( Table 4). Explanations of the grain-size data are given next to Equation (2).

Predictive Model
Nowadays the most commonly used models are the Witczak model [11,32], the Hirsch model [33,34] and models using artificial neural networks (ANN models) [35][36][37]. New predictive models are also constantly being developed [38][39][40]. A new predictive model for calculating the stiffness modulus for the four-point bending beam method (4PBB) was developed based on a form of the Witczak-El-Badawy model equation for the dynamic-modulus method (DM), commonly used in countries designing pavements according to MEPDG. The results of the Witczak-El-Badawy model are presented in [41]. The equation of this model is presented as Equation (1).

Predictive Model
Nowadays the most commonly used models are the Witczak model [11,32], the Hirsch model [33,34] and models using artificial neural networks (ANN models) [35][36][37]. New predictive models are also constantly being developed [38][39][40]. A new predictive model for calculating the stiffness modulus for the four-point bending beam method (4PBB) was developed based on a form of the Witczak-El-Badawy model equation for the dynamicmodulus method (DM), commonly used in countries designing pavements according to MEPDG. The results of the Witczak-El-Badawy model are presented in [41]. The equation of this model is presented as Equation (1). log E DM = 0.02 + 0.758· |G * | −0.0009 where: E DM -stiffness modulus determined by dynamic-modulus test [10 5 (1)) has undergone some modifications, also aimed at simplification. Some quotients (four) with the squares of the variables characterising the aggregate grain-size distribution were removed from the Witczak-El-Badawy model equation but other quotients with the variables characterizing the aggregate-size grain distribution were still left. The decision to remove some quotients from the equation was based on the results of the sensitivity analysis carried out for the Witczak-Bari model presented in [11]. This analysis shows that the variables with the smallest impact on the equation value are precisely the variables characterizing the aggregate grain size. In addition, the grain-size and air-voids content of HMA are the dependent variables and changes in grain size are partly taken into account by changing the air-voids content in HMA. Removing the quotients from the equation significantly simplifies both the model equation and the process of optimising the coefficients of the modified equation. Furthermore, the simplified model with 17 coefficients gave better fitting quality parameters than the Witczak-El-Badawy model equation with 21 coefficients. The equation of the modified model is presented as Equation (2). After these modifications, the model was named model A. The types of sieves used to characterize the grain size of HMA were also modified. In the original Witczak-El-Badawy model equation, sieves typical of AASHTO-associated countries, such as the sieve with a mesh side dimension of 0.075 mm, were used. This set of sieves has been replaced by a set used in some European countries (e.g., in Poland), with a sieve mesh side size as close as possible to the sieves being replaced.
where: As already stated, 471 values of the stiffness modulus E 4PB and the phase angle Φ representing the results of determinations made on the eight HMA mixtures were used to develop model A. The range of variability of the input data for the optimization process, i.e., the variables in the models, is given in Table 5. Optimization of the coefficients in the equations of the models was carried out using the lsqnonlin command in Matlab, which returns the solution of the nonlinear least-squares method task. The least-squares method is a method of approximating a function of a given type, to a set of empirical points. The method consists in choosing the parameters of the function being approximated in such a way that the sum of squares of the deviations of the empirical points from the value of this function is as small as possible.
Optimization was carried out both by specifying boundary conditions and without specifying them. In the case of model A, the solution depended on boundary conditions, similar to the authors of this model [32]. The optimization process for model A was carried out from a very large number of starting points determined in different ways, also by drawing the values of boundary conditions from a given interval. The evaluation of the solution took into account the value of the sum of least squares (∑e i 2 ), the limiting value of relative errors (P95(re)) and the distribution of absolute errors ( Table 6). As one can read in the work [32], it was discussed that in order to be considered a valid model, the distribution of absolute errors, after discarding outliers (nout), should be a normal distribution, and it was considered that the normality of the distribution should be confirmed by both histogram and statistical tests. The coefficients calculated as a result of the optimization procedure were presented in Equation (3)   Parameters used in Equation (3) were explained under Equation (2). Evaluating the model according to the graphs (Figures 16-18) and Table 7, it can be concluded that it has high precision and low bias (indicators explained in [32]). In Figure 16, it can be seen that for the HMA 7 and HMA 8 mixes, the determined E 4PB stiffness moduli are greater than the calculated ones, while the opposite is true for the HMA 2 and HMA 6 mixes. The influence of the type of aggregate and, more specifically, the silica content (SiO 2 ) in the aggregate, which affects the adhesion of the asphalt to the aggregate, can be seen here. For mixtures with more alkaline aggregates (HMA 7 and HMA 8), the determined stiffness moduli of E 4PB are larger than calculated. Therefore, at the stage of further work on the predictive models, the introduction of an additional variable characterizing the silica content of the aggregate used to make the HMA can be considered in the models. In Figure 17, where the test temperatures T are marked, it can be clearly seen that for high temperature, lower values of the stiffness modulus E 4PB were obtained, and for low temperature, higher values of the stiffness modulus E 4PB were obtained, confirming the generally known relationship and confirming the correctness of the developed model A. It can be noted, however, that for the analysed HMA, the values of stiffness moduli E 4PB obtained at 10 • C are only slightly smaller than those obtained at 0 • C. In Figure 18, where the loading frequencies f are marked, no significant relationship is noticed.
To further confirm the validity of the developed model for calculating the E 4PB stiffness modulus, a normality analysis of the distribution of the model's absolute errors was performed. The analysis of the normality of the distribution was performed by plotting histograms and performing multiple statistical tests for the normality of the distribution, with a significance level of α = 0.05. The D'Agostino-Pearson test was used as the binding test for the normality of the distribution due to the fact that the power of this test classifies at the medium level. Before analyzing normality, deviating absolute error values were rejected using the Hampel statistical test.
After applying the Hampel statistical test to eliminate outliers from the set of absolute errors for model A, 65 outliers were eliminated. The histogram ( Figure 19) indicates that the distribution of errors may be a normal distribution. The results of the applied statistical tests for the normality of the distribution are summarized in Table 8. The distribution of errors can be considered normal when the value of the p statistic > 0.05. According to all applied tests, the distribution of absolute errors can be considered a normal distribution. Almost all outliers, with few exceptions, are values calculated for results obtained at 0 • C.
Evaluating the model according to the graphs (Figures 16-18) and Table 7, it can be concluded that it has high precision and low bias (indicators explained in [32]). In Figure  16, it can be seen that for the HMA 7 and HMA 8 mixes, the determined E4PB stiffness moduli are greater than the calculated ones, while the opposite is true for the HMA 2 and HMA 6 mixes. The influence of the type of aggregate and, more specifically, the silica content (SiO2) in the aggregate, which affects the adhesion of the asphalt to the aggregate, can be seen here. For mixtures with more alkaline aggregates (HMA 7 and HMA 8), the determined stiffness moduli of E4PB are larger than calculated. Therefore, at the stage of further work on the predictive models, the introduction of an additional variable characterizing the silica content of the aggregate used to make the HMA can be considered in the models. In Figure 17, where the test temperatures T are marked, it can be clearly seen that for high temperature, lower values of the stiffness modulus E4PB were obtained, and for low temperature, higher values of the stiffness modulus E4PB were obtained, confirming the generally known relationship and confirming the correctness of the developed model A. It can be noted, however, that for the analysed HMA, the values of stiffness moduli E4PB obtained at 10 °C are only slightly smaller than those obtained at 0 °C. In Figure 18, where the loading frequencies f are marked, no significant relationship is noticed.          where the values are so small that even a small difference will result in a significant relative error. Medium relative error values were found in the results of tests carried out at 0 • C. The smallest relative error values were found at 10 • C and 20 • C. The relative error values at 10 • C are within approx. 30% and at 20 • C within approx. 20%, which can be considered satisfactory.

Concluding Remarks
According to the study and its analysis, further search for a better formulation of the model equation may not yield significant results. However, the introduction of an addi-

Concluding Remarks
According to the study and its analysis, further search for a better formulation of the model equation may not yield significant results. However, the introduction of an additional variable into the models to determine the silica content of the aggregate used to make the HMA or variables characterizing the shape of the aggregate grains could be considered. Further efforts to increase the accuracy of predictive models should focus on reducing the measurement uncertainty of the variables used in the models, which would be in line with the work described in [34].
It is also very important to assess predictive models in terms of relative errors. It can be concluded that in the literature on these models, information containing relative error analysis is very scarce. Among many publications analyzed, the authors found information on relative errors only in the article [34], where the figure includes simple ones specifying a relative error of re = 50%. Model A developed in this paper has just such a value for relative errors (P95(re) = 50%). The absolute errors of model A have a normal distribution similar to the model described in the work [11]. The correlation between the values obtained from the laboratory tests (Measured E 4PB ) and the values obtained from the model (Predicted E 4PB ) is very high, the coefficient of determination being R 2 = 0.936.
Model A shows the greatest inaccuracy in calculating the stiffness modulus E 4PB at 0 • C, confirming the results of the work of Bari and Witczak [32]. Similar to the mentioned researchers, the authors look for physical hardening effects here.
The authors are of the opinion that, in order to increase the accuracy of model A, a number of measures are needed to reduce the uncertainty of both the input variables of the model in question (especially a reduction in the uncertainty of the shear modulus |G*| of the bitumen) and the stiffness modulus E 4PB of HMA itself. A more precise procedure for determining the stiffness modulus E 4PB would contribute to a reduction in uncertainty. The conditions for thermostating the specimens prior to the low-temperature determination of the stiffness modulus E 4PB of HMA need to be defined more precisely in order to further reduce the possible influence of physical hardening. The algorithm for calculating the stiffness modulus E 4PB of HMA also needs to be clearly defined.
Using model A, it is proposed that, in addition to calculating the value of the HMA stiffness modulus E 4PB , it is also proposed to calculate the limiting values at the significance level α = 0.05 within which the result of the determination should fall, i.e., the uncertainty interval. These limits are proposed to be calculated using the P95(re) statistic. This would be done by subtracting and adding 50% of this value from the calculated value of the stiffness modulus E 4PB of the HMA; in this case, it would be known in which interval with 95% probability the laboratory-determined value of the stiffness modulus E 4PB should fall.