Computational generation of multiphase asphalt nanostructures using random fields

This study presents a novel methodology to generate computational replicates of nanostructures of multiphase materials, such as asphalt binders, by integrating image analysis techniques with stochastic random field (RF) modeling. Image analysis techniques are used to identify and segment nanostructure images obtained by atomic force microscopy, while RF is used to model the spatial distribution of their material properties. The results of this process are images showing probable arrangements of nanostructures with stochastic material properties that replicate the experimentally obtained images. The computationally generated nanostructures are then used as inputs in a finite element model to evaluate the effect of heterogeneity on their mechanical response. The efficacy of the developed approach is demonstrated through simulations of asphalt binders’ nanostructures, which reveal novel insights regarding their nanoscale mechanical behavior and response. The FE simulations provided the link between the distribution of nanoscale properties of asphalt binders and variations in their mechanical response. The application of this methodology expands the body of knowledge beyond the deterministic analysis of asphalt binders toward probabilistic analysis and uncertainty quantification that considers their heterogeneous, multiphase structures. Consequently, the methodology can be used to design multiphase materials, such as asphaltic blends, with tailored properties and enhanced performance.


INTRODUCTION
The variable and heterogeneous nature of multiphase materials impose uncertainty on their mechanical behavior as a consequence of variations in their constituent properties. An example of a multiphase material is asphalt binder, a material that consists of different nanoscopic phases that exhibit variable properties (Allen et al., 2012;Das et al., 2013). The intrinsic heterogeneity of such materials at the nanoscale is translated into variations in their mechanical response and when used as part of civil works, such as road structures, their field performance, and durability. On top of that, the widely adopted addition of asphalt modifiers introduces further heterogeneity to the material due to the presence of new constituents/phases Roja et al., 2022). It is essential to account for the effects of these variabilities to understand their behavior and predict performance; computational tools are especially well-suited for this task. Traditionally, the performance of asphalt binders is evaluated through rheological and mechanical tests at the bulk scale. However, nanoscale tests have gained popularity in recent years due to their capabilities in identifying the subcomponents of asphalt binders. One method that has been widely used to characterize asphalt binders at the nanoscale is atomic force microscopy (AFM). Several studies (Holleran et al., 2020;Roja et al., 2020;Wang & Liu, 2017) have used AFM to image asphalt binders and measure their mechanical properties at the nanoscale. AFM characterizes nanosurfaces. However, to maintain consistency with the general use in the literature, the term nanostructures will be adopted herein. Images produced by AFM portray the spatial arrangement of the nanostructures and their mechanical properties. Therefore, one can take advantage of the produced images to be used as key inputs in computer-aided mechanical numerical models to evaluate the response of asphalt binders. Such a method can help expand the knowledge of the performance of binders and relate their bulk performance to nanostructures and constituents, enhancing the quality of performance prediction. Indeed, with the aid of AFM, researchers have developed computational mechanical models simulating asphalt binder nanostructures. For example, researchers including (Du, Liu, Ganchev et al., 2021;Jahangir et al., 2015;Li-yan et al., 2020;Macedo et al., 2020) developed finite element (FE) models to study the internal stress distribution in binders' nanostructures subjected to external loads using geometric configurations and material properties obtained from AFM testing.
Generally, the development of mechanical FE models with only one geometric configuration (or nanostructure distribution) limits it to deterministic analysis. Moreover, assigning a single material property value to each phase of the FE model indicates the homogeneity of their material properties (Adeli & Kumar, 1995;Yu & Adeli, 1993). In reality, many materials, including asphalt binders, exhibit heterogeneity at the nanoscale (Aljarrah & Masad, 2020;Aljarrah et al., 2022;Holleran et al., 2020). This heterogeneity is demonstrated by the spatial variability of their nanostructures and material properties. Therefore, it is not realistic to assume deterministic and homogeneous conditions when modeling asphalt binders. For example, the mechanical response of such models may depend on the spatial arrangement of the nanostructure relative to the applied boundary conditions. Therefore, to properly characterize the mechanical response of the material, it is crucial to account for variations in its nanostructure distribution and the inherent heterogeneity of its multiphase properties.
In order to account for such spatial variability, a large number of representative AFM experiments should be carried out, which is an expensive and time-intensive process. These limitations hinder the feasibility of obtaining a large number of test results. One solution to overcome such limitation is by randomly generating probable computational replicates of the experimental results obtained by AFM. Such replicates should account for the spatial variability of the different phases to expand insights beyond the deterministic assumption of computational mechanical models. Many studies at the macro-and microscale levels have been performed to achieve a similar goal. In the area of asphalt mixtures, for example, (Aragão et al., 2017;Castillo et al., 2015;Klimczak & Cecot, 2020) generated computational replicates of coarse aggregate particles and air voids within asphalt mixtures to serve as inputs in computational mechanics models. Yet, no study has addressed the challenge of generating computational replicates of asphalt binders at the nanoscale to obtain probable representations of the materials' nanostructure, mirroring those obtained by experimental imaging techniques.
Likewise, stochastic modeling can be implemented to extend insights beyond the homogeneity of material properties. When used for these purposes, stochastic modeling assigns random yet spatially correlated values of material properties across a predefined space. One technique to achieve this goal is random field (RF) theory (Vanmarcke, 1983(Vanmarcke, , 2010. RF theory is a stochastic process that generates vectors of probable and spatially correlated random numbers within a spatial domain. The random numbers represent a variation of a specific parameter within the spatial field of interest. An example of such a parameter is the air void distribution within asphalt mixtures. Several studies have implemented RF to model heterogeneity in pavement materials Lea & Harvey, 2015;Manrique-Sanchez et al., 2020). Combining stochastic modeling with FE simulations offers the advantage of studying the effects of material properties on the mechanical response of materials without the need to build complex FE models. By doing so, the heterogeneity is captured while keeping the computational cost of the probabilistic model affordable. Moreover, such studies enable characterizing the expected or mean value of the mechanical response along with a measure of variability such as a standard deviation, which aids the probabilistic analysis and uncertainty quantification of the response of the materials .
Nonetheless, a number of previous studies on stochastic modeling of heterogeneous materials have attempted to represent multiphase materials using a set of statistical functions without differentiating among the phases. Effectively, this process generates a homogenized, single-phase that may not represent the complex response associated with the multiphase materials Hun et al., 2019;Stefanou et al., 2017). Other studies utilized microstructure generation software to produce structures with predefined geometric shapes Hou et al., 2016;Kim et al., 2016). Finally, some studies utilized non-Gaussian RF models or Gaussian RF models with nonstationarity processes and multiple-point correlation functions (Fuglstad et al., 2015;Montoya-Noguera et al., 2019;Xu & Gardoni, 2020). Such models are computationally complex when it comes to determining the parameters of functions that describe the distribution of multiphase materials. In light of these limitations, our method is proposed in this context to try and address them and complement a gap in the state of the art.
This study proposes a novel methodology based on integrating image analysis techniques with stochastic RF modeling to generate two-dimensional (2D) random nanostructures of multiphase materials exhibiting heterogeneity or spatially variable properties. Such variability, which is induced by the generated nanostructures, can be translated into variations in materials' mechanical response, thus providing insights regarding the uncertainty of their behavior in the field (i.e., reliability of the structure). Asphalt binders were used as an example of such materials. The randomly generated nanostructures were used as inputs in computer-aided mechanical models to evaluate the role of their heterogeneity on the mechanical response of asphalt binders. Computational mechanical modeling helps in understanding asphalts' local stress distribution as a function of the spatial arrangement of their components and local properties (Jahangir et al., 2015). The methodology also enables detailed analysis beyond deterministic analysis of FE modeling and the assumption of homogenous material properties. The method also aids in better designs of multiphase materials, such as asphaltic blends with tailored properties and enhanced performance. The added value of integrating RF with computational mechanics modeling is the possibility of assigning unique values of material properties to each element in the model mesh, which helps to account for the effects of material heterogeneity (Castillo & Al-Qadi, 2018). In addition, creating this virtual testing environment promotes sustainability by conserving the resources required to carry out laboratory tests and enables the evaluation of hypothetical scenarios that would be impractical in a laboratory environment.

RF THEORY
RF theory is a branch of the general stochastic theory (Vanmarcke, 1983), which is a process that generates a set of arrays of probable and spatially correlated random numbers within a spatial domain (Adler & Taylor, 2007;Vanmarcke, 2010). The random numbers represent variations in a specific parameter within the spatial field. Spatial correlation means that neighboring locations within the field are more likely to have closer values of the parameter of interest. One generation of an RF (i.e., a geometry of spatially correlated random numbers of the selected parameter within the predefined space) is called a "realization." In each realization, the generated arrays contain information regarding the magnitude, coordinates, and correlation characteristics of the parameter of interest. The covariance matrix decomposition (CMD) technique (El-Kadi & Williams, 2000) has been proven to be an efficient RF method to generate normally distributed 2D RFs Narsilio, 2006). The set of input parameters of this technique includes calculating the mean (μ) and standard deviation (σ) of the parameter of interest, estimating the field correlation lengths (L x and L y in a Cartesian x-y coordinate system), determining the dimensions of the field (n × m), and choosing a suitable correlation function (f) to regulate how values are correlated within the field. Herein, the set of input values (i.e., μ, σ, L x , L y , n, and m) were determined by analyzing the distribution of the parameter of interest. Regarding f, an exponential correlation function is adopted based on the work of . The horizontal and vertical correlation lengths (L x and L y ) determine how values are correlated within a spatial range in each direction; values within this range are expected to be highly correlated; within the range of a correlation length meaning that the likelihood of finding similar values increases. The "expeditive" method (Stuedlein et al., 2012;Vanmarcke, 1977) was selected in this study to estimate their values. The method suggests that the correlation length can be estimated as eight-tenths of the average distances between the intersections between the profiles of the parameter values and their linear trend.
To produce one random realization, the values of the RF parameter Y can be calculated using the following "Equation (1)": where Y is a realization of the RF, which comprises a vector of the parameter of interest. A is the autocorrelation matrix, which will be elaborated shortly. G is a vector of standard normally distributed random values, that is, G∼N[0,1], and μ is the mean vector of the parameter of interest. The autocorrelation matrix, A, is calculated by applying the Cholesky decomposition: where C is the covariance matrix of the field. Cholesky decomposition is a practical technique for decomposing Hermitian positive definite matrices into the product of a lower triangular matrix and its conjugate transpose. The covariance matrix C comprises information regarding the parameter standard deviation, along with its governing spatial correlation function, expressed as where f i,j is the correlation function calculated for each pair of points (i,j) in the spatial domain. An exponential correlation function was adopted in this work to correlate the values of the parameter of interest, expressed as follows: where , and , are the horizontal and vertical distances between the pair of points i and j within the space of interest. In this study, the implementation and generation of RF realizations are conducted in the MATLAB environment.

Characterization of nanostructures using AFM
In order to generate computational replicates of asphalt binders' nanostructures, experimental procedures are initially carried out to characterize the nanostructures by imaging their arrangement and measuring their nanomechanical properties. To this end, the peakforce quantitative F I G U R E 1 Typical atomic force microscopy (AFM) image obtained by the peakforce quantitative nanomechanical mapping test for a modified asphalt binder. The color bar represents the variation in the Derjaguin-Muller-Toporov (DMT) modulus across the sample surface nanomechanical mapping (PFQNM) test was carried out using Bruker's Dimension Icon AFM equipped with a Nanoscope. The PFQNM test produces images showing the spatial arrangement of samples' nanostructures while measuring their Derjaguin-Muller-Toporov (DMT) elastic modulus (Derjaguin et al., 1994). The DMT modulus, derived from the DMT contact mechanic model, is a modified version of the Hertzian model with the adhesive forces initiating between the two bodies. The DMT model is given by the following "Equation (5)": where DMT is the modulus value, which corresponds to Y in "Equation (1)," P is the force applied to the sample surface, a is the adhesion force between the AFM tip and the sample surface, r is the tip radius, and Δ is the deformation of the sample surface. In this study, the PFQNM test used the RTESPA-300-30 tip with a radius of 30 nm and a spring constant of 40 N/m. The applied peak force was set to 50 nN, and images were produced with dimensions of 10 × 10 μm and a resolution of 2000 × 2000 pixels. Figure 1 shows a typical AFM image obtained by the PFQNM test for asphalt binder modified with 5% low-density polyethylene (LDPE). This image, which was selected as an example throughout this section to explain the proposed methodology, shows the spatial arrangement of the nanostructures and their DMT modulus across the sample surface. The lower-modulus circular nanostructures will be referred to as "dispersions," or the dispersed phase, and F I G U R E 2 Computationally generated nanostructures assuming homogenized phases the background medium containing the dispersions will be referred to as the continuous phase. The color bar indicates the range of the DMT modulus distribution across the sample surface. It should be noted that the mean and standard deviation values of the DMT modulus distribution in Figure 1 are 1.145 and 0.069 GPa, respectively. As shown in the following sections, RF modeling presented in Section 2 is used to represent the spatial variation of DMT moduli of asphalt binders. More information about the binders characterized in this study is presented later in this paper.

Computational generation of single-phase nanostructures of DMT modulus by means of RF
As an initial attempt, the proposed methodology intended to model the experimentally captured images using a single RF to represent the spatial distribution of the modulus of the entire image. Here, values of μ = 1.145 GPa and σ = 0.069 GPa were used based on the experimental image to generate an RF realization for the image as a whole. By doing so, all nanostructures are homogenized into one single phase. L x and L y values were iteratively chosen as 7000. The generated image of the nanostructures is shown in Figure 2. The RF has dimensions of 100 × 100 pixels, corresponding to a pixel size of 100 nm. The resulting mean and standard deviation values are 1.178 and 0.052 GPa, respectively. Two key observations can be drawn from Figure 2. One, the distinct dispersions shown in Figure 1 are not as pronounced in the computational replicate. They look more "diffused" in the continuous phase. Second, the range of the color bar significantly differs from the original image. Although both images have proximate average and standard deviation values of the DMT modulus, they represent different structures, as the variability of the original structure is not properly modeled. These differences are due to the limitations of the single-phase assumption. Thus, modeling the image using single-phase nanostructures does not represent the original nanostructure distribution and properties, and a different approach should be adopted.
As mentioned earlier, the proposed methodology applies the CMD technique to generate spatially correlated RFs following a Gaussian distribution. Moreover, the method relies on an exponential two-point correlation function and assumes stationary correlation lengths throughout the field. Nonetheless, the spatial distribution of the DMT modulus shown in Figure 1 does not follow a Gaussian distribution, and the presence of dispersions throughout the continuous phase will yield nonstationary correlation lengths. This would adversely affect the accuracy of the RF model and produce unrealistic replicates of the nanostructures as noted in Figure 2. To overcome such issues, one can adopt higher-order correlation functions based on nonstationary multiple-point correlation lengths. In addition, non-Gaussian RFs can be used to generate the values of the parameter that will be assigned to the different locations of the field (Montoya-Noguera et al., 2019). However, both approaches will impose higher computational expenses. Therefore, an alternative solution is proposed in this study by integrating the capabilities of image analysis techniques with Gaussian RF models; this approach significantly increases modeling accuracy and simplifies the computational processes.

Computational generation of multiphase nanostructures by means of RF
A novel methodology based on image analysis techniques and stochastic RF modeling is proposed in this study to replicate the nanostructures of multiphase materials. An illustrative flowchart explaining the methodology is shown in Figure 3. In summary, experimentally obtained images are analyzed to extract the value of the DMT modulus at each pixel and identify nanostructures and then segment them into phases. Next, statistical analyses of the distribution of the DMT modulus are conducted for each phase in order to calculate the set of inputs required to build RF models. Subsequently, RF realizations are generated to model the spatial distribution of the DMT modulus for each phase. Finally, the generated phases are randomly arranged in one spatial domain, creating a computational replicate of the experimental image. The process of generating RF realizations for individual phases and then combining them is repeated to generate many replicates of the F I G U R E 3 Flowchart of the proposed methodology to produce random nanostructures of multiphase materials experimental image. The AFM image of a modified asphalt binder presented in Figure 1 was selected to illustrate the steps of the proposed method. However, a deeper analysis of the application of the method to evaluate the impact of the nanostructure heterogeneity on the mechanical response of asphalt binders is presented in a later section. A detailed description of each step is explained next.

Image analysis
The first step in the proposed methodology entails analyzing the AFM image to extract the value of the DMT modulus at each pixel. This is achieved by converting the image to grayscale and then calibrating the grayscale values based on the color bar shown in Figure 1. The grayscale image is shown in Figure 4. Next, image analysis techniques are implemented to identify nanostructures and segment them into phases of similar features (e.g., geometric shape or color). To do so, the Image Segmenter App in the MATLAB environment was utilized. Appropriate image segmentation tools are employed, such as automatic detection tools (Adaptive Thresholding, Auto F I G U R E 4 AFM image after converting to grayscale Cluster, and Find Circles), manual detection tools (Flood Fill), and refinement masks (Fill Holes, Dilation, Erosion, Opening, and Closing). The segmentation process is carried out iteratively since no ideal segmentation tool would accurately identify all types of features. Figure 5 shows a binary mask of the AFM image after applying image segmentation tools. The dispersions are F I G U R E 5 AFM image after applying image segmentation tools shown in white color, while the continuous phase is represented by the black background. The dispersions are features or clusters that have specific geometric shapes and color shades, where color represents the DMT modulus value. By segmenting the nanostructures into distinguishable phases and obtaining modulus values across each phase, the set of inputs required to construct the RF models can be calculated.

Statistical analysis of image phases
At this stage, the inputs required to build the RF models are calculated. This is done by constructing the frequency distribution of the DMT modulus for each phase. As an example, Figure 6 shows the frequency distribution of the DMT modulus of the largest dispersion in Figure 1. From the n × m matrix containing modulus coordinates and values, μ and σ were calculated. Likewise, the expeditive method (Stuedlein et al., 2012) is applied to estimate the values of L x and L y . As mentioned earlier, the CMD technique can generate normally distributed RFs; however, if the distribution of the parameter of interest exhibits skewness, a lognormal transformation, that is, k = log x Y can be applied. In this study, a lognormal transformation of base two was applied before calculating RF inputs. The value of two was chosen iteratively. The frequency distribution shown in Figure 6 was constructed after applying the transformation.

RF modeling of the spatial distribution of the parameter of interest
The CMD technique can now be applied to create RF realizations to model the spatial distribution of the DMT F I G U R E 6 DMT modulus frequency distribution of the largest dispersion previously shown in Figure 1. The left inset shows the original dispersion after being segmented. The inset to the right shows the values of the DMT modulus within this dispersion modulus for each phase. This is done by substituting the calculated values of μ, σ, L x , L y , n, and m into "Equations (1) through (4)." Initially, the horizontal and vertical distances between each pair of pixels (θ i,j x and θ i,j y ) are substituted in the correlation function shown in "Equation (4)" along with the values of L x and L y calculated in the previous step. Then, the covariance matrix (C) is constructed using "Equation (3)." Subsequently, the Cholesky decomposition is applied to matrix C as shown in "Equation (2)," to obtain the autocorrelation matrix (A). Afterward, matrix A is multiplied by a vector of standard normally distributed random values (G) and scaled by the mean vector of the DMT modulus (μ), producing one realization of the RF (Y) for that phase. The methodology starts by generating an RF realization to model the spatial distribution of the DMT modulus across the continuous phase to serve as the background in the final computational image. Next, there are two options to model the distribution of the modulus across the dispersions: The first option is to generate RF realizations for each dispersion. The second one is to use original DMT modulus values as measured by AFM. The second option offers the ability to evaluate the effect of randomness in nanostructures' spatial arrangement without varying their modulus values.

Generation of computational images
Finally, the generated dispersions are randomly arranged across the generated continuous phase following a uniform distribution, thus creating a computational replicate of the experimental image obtained by AFM. The process F I G U R E 7 Computationally generated images of resolution 100 × 100 pixels, replicating the experimental AFM image shown in Figure 1 of generating RF realizations for individual phases (Step 3.3.3) and then randomly placing them in one spatial domain (Step 3.3.4) is repeated to generate different replicates of the experimental image. Figure 7 shows nine computational replicates of the AFM image shown in Figure 1. These images exhibit probable arrangements of the nanostructures and spatially correlated DMT modulus values. The resolution of the generated images is 100 × 100 pixels. The entire code, which comprises Steps 3.3.1 through 3.3.4, was executed on a workstation with 12 GB of available memory. While the execution time for the first three steps was around 530 s, the image synthesis step typically took less than a second.

Validation of the proposed methodology
To validate the accuracy and repeatability of the proposed methodology, frequency distributions of the DMT modulus were compared for (1) the experimentally imaged sample (Section 3.1), (2) a computational replicate assuming single-phase nanostructures, (Section 3.2), (3) a computational replicate assuming multiphase nanostructures without a lognormal transformation, (Section 3.3), and (4) 25 computational replicates assuming multiphase nanostructures with a lognormal transformation (Section 3.3). More than 25 replicates can be easily added. Table 1  table indicates that the mean of the computational images matches the original mean within one standard deviation. However, the closest mean value was for the average of the 25 replicates, indicating the accuracy of this method. A sole comparison of the means and standard deviations of the modulus is insufficient. Other features, such as the shape of the distributions and the overall appearance of the fields, must also be representative. Figure 8 shows the frequency distributions of the DMT modulus for the four cases. Clearly, the assumption of single-phase nanostructures produced the least accurate distribution. Moreover, the figure indicates that applying the lognormal transformation greatly enhanced modeling accuracy, especially at the tails of the distribution. Finally, by comparing the distributions of the 25 computational replicates (shown in gray solid lines), the repeatability of the proposed methodology is confirmed by the proximate distributions.
To further assess the statistical similarity between the experimental and computational replicates, the distributions in Figure 8 are also presented as boxplots in Figure 9. As noted, all computational replicates showed F I G U R E 9 Boxplots of the distribution of the DMT modulus for the experimental image and computational replicates with different assumptions overlapping interquartile ranges with the experimental one. Nonetheless, the highest similarity was spotted for the computational replicates obtained by the proposed methodology (multiphase-log). The proposed methodology managed to replicate the experimental data even beyond the maximum and minimum values. This reinforces the accuracy of the method. On the other hand, Figure 8 exhibits an interesting observation, where narrower distributions with lower peak heights were observed for the 25 replicates, compared to the experimental one, which aligns with the shorter range in Figure 9. This can be attributed to a large number of outliers in the experimental image, especially above the maximum value. These outliers could not be accurately modeled. Nonetheless, the proposed methodology captures the central portion of the distributions efficiently. However, further improvements are required to detect the full variability of the distributions.

Materials and sample preparation
To demonstrate the effectiveness of the proposed methodology, the effect of the spatial variability of the nanoscale properties on the mechanical response of asphalt binders was evaluated. Nanoscale properties refer to the arrangement of nanostructures and the spatial distribution of their DMT modulus. To this end, three different samples were prepared. An unmodified binder classified as Performance Grade 64-22 (PG64-22) was used as a control sample and used to prepare the other two modified samples. These other samples contained LDPE at 3% and 5% dosages by weight of the binder and are referred to as 3 and 5L, respectively. Both samples 3 and 5L had a high-temperature performance grade of PG70s. There is growing interest in the use of plastics as asphalt modifiers to achieve sustainability goals and because some types of plastics, such as LDPE, have been shown to enhance the high-temperature performance of asphalts Nuñez et al., 2014). The LDPE was added to the base binder and blended using a high shear mixer at a temperature of 185 • C for an hour. Subsequently, for AFM specimen preparation, the heat and cast method was adopted (Aljarrah & Masad, 2020), where approximately 150 mg of the hot blend was smeared on a glass slide and reheated for another 3 min to ensure a smooth surface. The samples were then cooled to ambient temperature (∼25 • C) and contained for 1 day before carrying out AFM testing.

Generating computational replicates of AFM images using the proposed methodology
The PFQNM test was carried out by means of AFM to image the samples' nanostructures and measure their DMT modulus values. The captured images are shown in Figure 10. As noted, the bee-like structures (dispersions) in sample PG64-22 are in fact ripples of peaks showing the highest modulus values and valleys showing the lowest. The continuous phase exhibits values varying between the two extremes. In samples 3 and 5L, the dispersions exhibit lower modulus values compared to the continuous phase, but in both samples, the average modulus values were higher than the PG64-22 sample, indicating higher stiffness.
After imaging the prepared samples by AFM, images were converted to grayscale and calibrated according to the color bars shown in Figure 10. Subsequently, images were segmented into distinguishable phases. Figure 11 shows AFM images after segmentation, and the dispersions are shown in white. The dispersions in sample PG64-22 are bee-like structures, and those in samples 3 and 5L are circular clusters.
Next, the inputs required to create the RF models were calculated. A lognormal transformation of base two was applied to the distribution of the DMT modulus before calculating μ, σ, L x , and L y . Subsequently, the spatial distributions of the DMT modulus across the continuous phases in all samples were modeled by RFs. Unlike samples 3 and 5L, the dispersions in sample PG64-22 were not modeled by RFs. Rather, the original modulus values were assigned. Finally, the generated dispersions were randomly placed across the generated continuous phase, producing computational replicates of AFM images previously shown in Figure 10. This process was repeated 25 times to ensure repeatability. Figure 12 shows some of these computational replicates. The resolution of these images is 100 × 100 pixels.

FE simulation of the mechanical response of nanostructures
The FE method is a tool used to acquire numerical solutions to problems that are tedious or costly to solve by analytical techniques. Therefore, FE models simulating the experimentally imaged and computationally generated nanostructures were developed to evaluate the role of their heterogeneity on the mechanical response of asphalt binders. This allows quantifying not only the expected value (mean) of their mechanical response but also offers a measure of its variability (e.g., standard deviation) to the analysis. Having information about the variability of the mechanical response of asphalt binders provides insights regarding the uncertainty of their behavior in the field (i.e., reliability of the structure). To this end, the commercial FE software ABAQUS was used to develop computational mechanics models simulating the experimental and computational nanostructures. The dimensions of the developed models are 10 × 10 μm, which were divided into 100 × 100 divisions, similar to the resolution of the images. The model was meshed into 100 × 100 elements as well, allocating one element per division. In other words, the size of the FE mesh is the same as the resolution of the images. The mesh size was selected after conducting a sensitivity analysis of its effect on the model outputs as will be explained in the following section.
The element type was assigned as a 2D four-node plane stress shell element to reflect the thin thickness of AFM samples. The material property input, which is the elastic modulus, was assigned for each division/element in the model mesh based on the distribution of the DMT modulus. In other words, each element is assigned a unique DMT modulus value. In addition, Poisson's ratio was set to a value of 0.4 based on the recommendations of previous studies (Benedetto et al., 2007;Gudmarsson et al., 2015). The FE model simulates a strain-controlled uniaxial tensile test, where a nondeformable rigid body is attached to the top of the model while constraining the translation of the other side in all directions. A 1 μm displacement (10% tensile strain) was applied. Figure 13 shows an illustration of the FE model experiment for one replicate of each sample obtained by RF. The distribution of the internal stresses and strains can be acquired from such simulations.

Results and discussion
Sensitivity analysis was conducted to determine the effect of the FE model mesh size on the model response. To this end, different mesh element sizes, including 100, 50, 33, and 25 nm, were developed and analyzed. These dimensions correspond to allocating 1, 4, 9, and 16 mesh elements per individual division. The outcomes of the models are shown in Figure 14. The figure clearly shows that the Von Mises stress (VMS) distributions resulting from different mesh sizes are almost identical, especially at the central part of the distribution. Yet, slightly wider ranges are observed for smaller mesh sizes as indicated by the boxplots. Therefore, a value of 100 nm element size, corresponding to 1 mesh element per division, was selected to provide the highest computational efficiency required to simulate a large set of samples. Figure 15 shows the spatial distribution of von Mises stress across three computationally generated images of each type of binder. This analysis was conducted after eliminating the boundary effects by removing the stress values in elements within a 1 μm distance from the top and bottom of the model. Generally, the modified samples induced higher stress distributions, compared to the PG64-22 sample, and sample 5L showed the highest stress values. Figure 15 also reveals an interesting observation: Sample PG64-22 indicates that the highest stresses are located near the peaks of the bee-like structures, while the lowest stresses are located near their valleys. Such observation indicates that stress localization occurs at the bee-like structures due to the large variation in stresses, which may eventually lead to crack initiation.
The phenomenon of localized stresses in asphalt binders due to external loads was also observed in a study by (Jahangir et al., 2015). However, their FE model showed that stresses are localized at the continuous phase rather than the bee-like structures. They also showed that the lowest stresses emerge in the bee-like structures and the highest at the continuous phase, contradictory to the findings in this study. However, their FE model assumed homogenous phases, that is, one modulus value for the bee structures (peaks and valleys) and another modulus value for the continuous phase. In the current work, unique modulus values were assigned to each element in the model mesh based on the spatial distribution of the DMT modulus (elementwise assignment). In this way, modulus values at the peaks of the bee structures will be different than those assigned at the valleys. This shows the advantage of the proposed methodology in delivering novel insights regarding the behavior of nanostructures by introducing an accurate representation of variation in material properties, thus extending the insights beyond the deterministic and homogenous modeling conditions.
For samples 3 and 5L, it is shown in Figure 15 that stresses are more evenly distributed across the sample surface. This indicates that modification of these samples helped to develop a more uniform stress distribution and reduce stress localizations across the sample. Slight localization was observed near the dispersions due to differences in DMT modulus values between the continuous F I G U R E 1 3 Illustration of the finite element model experiment for one replicate of each sample obtained by random field F I G U R E 1 4 Frequency distribution and boxplots of von Mises stress for sample 5L for different mesh sizes phase and the dispersions. Such locations can be the weakest points to initiate cracking. The coefficient of variation (CV) for the VMS distribution was calculated for 25 computationally generated replicates for each binder. The CV values ranged between 8.1% and 9.1% for PG64-22, 4.7% and 4.9% for 3L, and 4.9% and 5.2% for 5L. Moreover, CV values ranged between 14% and 22% within the bee structures in sample PG64-22. This further proves that modification reduced the variability (increased uniformity) of the stress distribution in the modified binders. Figure 16 shows the von Mises stress distribution and elastic modulus distribution for the experimentally obtained nanostructures and 25 computationally generated samples, respectively. The overall trend aligns with findings from a previous work conducted by the authors of this study , in which they measured the bulk dynamic modulus of the same binders using a dynamic shear rheometer. This previous study showed that sample 5L produced the highest stiffness, F I G U R E 1 5 Spatial distribution of von Mises stress computationally generated samples with a resolution of 100 × 100 pixels, (a) PG64-22, (b) 3L, and (c) 5L F I G U R E 1 6 Von Mises stress distribution and elastic modulus distribution of the experimentally obtained (solid lines) and computationally generated (dashed lines) samples followed by sample 3L, and then PG64-22. Further, the proposed methodology succeeded in capturing the effect of the heterogeneity of nanostructures on the mechanical response of asphalt binders. In other words, the variability was translated into variations in the mechanical response distributions indicated by the dashed lines in both subfigures. Such variability can be quantified by adding a variability measure (such as standard deviation) to the analysis. In this way, the mechanical response is thus not defined by a single deterministic value but rather by a distribution or a confidence interval considering the effect of the uncertainty induced by material heterogeneity. These results clearly indicate the benefit of studying the response and variability in material properties at the nanoscale in order to understand material bulk behavior. Moreover, the distributions in Figure 16 also align with the observations in Figure 15, where flatter distributions with lower peak heights were observed in the modified samples. This further indicates that stresses are more evenly distributed across the surface of the modified binders, which is more favorable to reduce cracking potential.
Finally, it should be noted that the analysis of this application example utilized material properties based on the DMT modulus distribution obtained by the PFQNM test. The DMT modulus is an elastic modulus and therefore does not account for the viscoelastic nature of asphalt binders (Dong et al., 2020;. Future work should include such viscoelastic properties, which can be obtained utilizing the nanoscale dynamic mechanical analysis (nDMA) test using AFM (Aljarrah & Masad, 2020).

SUMMARY AND CONCLUSION
In this study, a novel methodology based on an image analysis technique and stochastic RF modeling is proposed to generate computational replicates of multiphase materials' nanostructures exhibiting heterogeneity or spatially variable properties. Asphalt binders were used to demonstrate the effectiveness of the methodology, where AFM testing was carried out to characterize their nanoscale properties, image analysis techniques were used to segment nanostructures into phases, and RF modeling was utilized to model the spatial distribution of their moduli. The computationally generated nanostructures were used as inputs in FE models to evaluate the effect of their heterogeneity on the mechanical response of asphalt binders. The following conclusions can be drawn: 1. The proposed methodology succeeded not only in generating representative computational replicates of multiphase material nanostructures but also in modeling the variability of their spatial arrangement and material properties. Additionally, the methodology is applicable to multiphase materials that exhibit different features and at different scales. 2. The variability in the nanoscale properties of multiphase materials is reflected in their mechanical behavior as proven by FE simulations. Such variability can be translated into variations in their mechanical response, thus providing insights regarding the uncertainty of their behavior in the field (i.e., reliability of the structure).
3. By applying the methodology to asphalt binders, interesting observations were revealed by probabilistically analyzing the results of the FE models. The analysis showed that the variation in stresses between the peaks and valleys of the bee-like structures in asphalt binders leads to stress localization, which could cause cracks to initiate and nucleate from these structures. On the other hand, a wider stress distribution in modified binders can be interpreted as an enhanced and more even stress distribution across the sample surface, reducing the chance of stress localization. Moreover, studying variations in the mechanical behavior of materials at the nanoscale can be beneficial in understanding their overall bulk performance. 4. This work demonstrated that accounting for the variability in inherently heterogeneous materials can lead to interesting and novel findings by expanding beyond deterministic analysis and homogenous conditions, adding value to experimental data. Moreover, this work promotes sustainability by conserving resources required to carry out expensive laboratory tests, as well as enabling the evaluation of hypothetical scenarios that would be impractical in a laboratory environment. 5. This study highlights the significance of probabilistic analysis and uncertainty quantification of asphalts' mechanical response. Nonetheless, the methodology implemented in this study is for 2D images. Future work should consider three-dimensional modeling of asphalts' nanostructures using reconstruction techniques to understand the bulk behavior of the material. Moreover, the developed FE models utilized the elastic properties of asphalts. Future research may include viscoelastic properties of such materials, which can be measured using nDMA, along with continuum damage modeling. This analysis would contribute to capturing the effects on nanostructure distribution and properties on binder resistance to damage. Finally, further refinements of the image processing steps and including more AFM images may improve objects detection, reduce subjectivity in the segmentation process, and enhance the representation of binder nanostructure distributions.

A C K N O W L E D G M E N T S
This work was supported by the Qatar National Research Fund (QNRF): NPRP11S-1128-170041. All the statements are those of the authors. Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing. Open Access funding provided by the Qatar National Library.