3D Microstructure Simulation of Reactive Aggregate in Concrete from 2D Images as the Basis for ASR Simulation

The microstructure of alkali-reactive aggregates, especially the spatial distribution of the pore and reactive silica phase, plays a significant role in the process of the alkali silica reaction (ASR) in concrete, as it determines not only the reaction front of ASR but also the localization of the produced expansive product from where the cracking begins. However, the microstructure of the aggregate was either simplified or neglected in the current ASR simulation models. Due to the various particle sizes and heterogeneous distribution of the reactive silica in the aggregate, it is difficult to obtain a representative microstructure at a desired voxel size by using non-destructive computed tomography (CT) or focused ion beam milling combined with scanning electron microscopy (FIB-SEM). In order to fill this gap, this paper proposed a model that simulates the microstructures of the alkali-reactive aggregate based on 2D images. Five representative 3D microstructures with different pore and quartz fractions were simulated from SEM images. The simulated fraction, scattering density, as well as the autocorrelation function (ACF) of pore and quartz agreed well with the original ones. A 40×40×40 mm3 concrete cube with irregular coarse aggregates was then simulated with the aggregate assembled by the five representative microstructures. The average pore (at microscale μm) and quartz fractions of the cube matched well with the X-ray diffraction (XRD) and Mercury intrusion porosimetry (MIP) results. The simulated microstructures can be used as a basis for simulation of the chemical reaction of ASR at a microscale.


Introduction
The cracking caused by internal chemical reactions, such as delayed etrringite formation and alkali silica reaction (ASR), is one of the major durability problems in concrete structures [1][2][3]. ASR has been viewed as a 'cancer' in concrete structures since it was firstly discovered and named by Stanton in the 1940s [4]. ASR starts from the dissolution of the reactive silica in the aggregate under the attack of hydroxide ions that come from the pore solution in the cement paste [5]. The dissolved silica ions either react with the diffused alkali ions inside the aggregate or diffuse into the cement paste and then react with the alkali ions. Where the products will be formed depends on a lot of factors. According to the research of [6][7][8][9][10], the mineralogical nature of the aggregate greatly affects the cracking pattern of the ASR-affected concrete. According to the study of Ponce and Batic [6], for rapid-reacting aggregates containing amorphous silica such as opal and vitreous volcanic rocks, ASR product formation and expansion occur at the paste-aggregate interface transition zone (ITZ). However, for slow-reacting aggregates containing defected crystal quartz, such as siliceous limestone, gel occurrence and expansion come mainly from the inside of the aggregate. In other words, the microstructure of the reactive aggregate partially determines the initial expansion sites caused by ASR in the concrete.
Developing predictive models is important since it is able to give suggestions about ASR prevention for the to-be-built structure or maintenance measures for the ASR-damaged structures. In the present predictive mechanical models, several models have considered the spatial reaction and initial expansion sites. The models of [11,12] simulated the mechanical degradation of the concrete induced by ASR. The reaction and expansion sites were represented by several expansive points or gel pockets distributed randomly in the aggregate or in both the aggregate and the ITZ. A recent mesoscopic model [13] has explored the influence of the initial expansion sites on the ASR cracking patterns. The expansion sites are assumed to be a reaction rim on the surface of the aggregate or to be some gel pockets distributed randomly inside sphere aggregates. These models are able to provide acceptable predictions about the cracking propagation or the structural performance at a mesoscale (mm) or macroscale (m), even though the expansion sites are simplified. However, the chemical reaction mechanism, such as when and where ASR products begin to grow, how the expansion sites evolve, the influence of the mineralogy of the aggregate on the expansion sites, and so forth cannot be captured or explored by these models. In other words, the real microstructure of the aggregate is needed if one tries to simulate ASR as realistically as possible or if one tries to build a model that can drive the physical response from the chemical reaction.
Our project is to build such a realistic model. The first problem is how to obtain the 3D microstructure of the aggregate. The most regular computed tomography (CT) or focused ion beam milling combined with scanning electron microscopy (FIB-SEM) is a direct method to obtain the 3D microstructure of a porous material [14]. However, unlike cement paste, whose representative elementary volume (REV) is around 100 × 100 × 100 µm 3 , the REV of a coarse aggregate is usually at a mesoscale (mm) [15]. Such a REV can be imaged with CT; however, it is presently challenging to resolve the pores and silica particles with sizes below 1 µm with this technique. With FIB-SEM, the required spatial voxel size can be met; however, the sample size of several micrometers is too small to be representative of the heterogeneous aggregates. Additionally, applying both of these methods to characterize the aggregates can be done [16]. Yet, such a workflow is time-consuming and expensive, as specific equipment with limited availability is required.
Alternatively, the 3D microstructure can be derived by simulating it based on the information from 2D images. This approach could deliver 3D information in a fast and inexpensive manner. A number of simulation approaches have been developed in the past [17], and among them, the autocorrelation function method and simulated annealing method [18] were extensively examined. The autocorrelation function method has been used and extended by many researchers [19][20][21] to simulate the 3D representative microstructure of materials. This method firstly generates Gaussian population noises in a domain and then filters the noises based on the 2D autocorrelation function to match the phase volume fractions and surface-area fractions on the 2D images.
Based on their experiences [20,22], the Gaussian filtering method and image analysis are utilized in this paper to simulate the 3D microstructure (mainly the spatial distribution of pore and reactive silica) of an alkali-reactive siliceous limestone from Belgium from 2D images. Instead of finding a REV for the aggregate at a mesoscale, five separated representative microstructures (each one with a size of 100 × 100 × 100 µm 3 ) with different pore fraction and reactive silica fraction following a true distribution in 2D images were simulated and upscaled to compose the microstructure of an irregular coarse aggregate in a simulated 40 × 40 × 40 mm 3 concrete cube. The simulated microstructure of the limestone is compared visually (phase distribution characteristics), numerically (ACF, particle size distribution), and experimentally (porosity, silica fraction) with the original 2D images, and X-ray diffraction (XRD) and mercury intrusion porosimetry (MIP) results.

Research Significance
A better understanding of the influence of the aggregate microstructure on ASR is essential to develop the materials' science-based design procedure for the newly-built concrete to prevent ASR. This study is an effort in that direction, where a 3D aggregate microstructure is simulated from 2D images providing an advantage to overcome the resolution limitation of a CT scan or FIB-SEM. The proposed method can be used by the ASR modelers to consider the realistic microstructure of the aggregate instead of assuming the aggregate to be a homogeneous sphere so that the fundamental mechanism can be simulated more realistically. Moreover, the method can be extended to other porous materials rather than aggregates. For example, the microstructure of cement paste can be reconstructed directly based on one representative 2D image.

Material and Experimental Methods
An alkali-reactive siliceous limestone from Belgium, which has been reported to cause damage to several concrete structures, such as bridges and dams [23], is used as the simulation object in the paper. Physical properties were determined firstly to distinguish the aggregate.
A 3D microstructure simulation is mainly based on the utilization of 2D quantitative information which can only be calculated from 2D microstructure images with a clear outline between particles. Based on this principle, optical petrography micrographs were firstly taken from thin sections of the limestone to identify its mineral compositions and their particle size. According to the thin section results, most of the minerals have a cryptocrystalline size (below 4 µm), which cannot be segmented by the Leica DMLP petrographic microscope (Leica, Wetzlar, Germany) with a maximum resolution of 200 mm in the laboratory. Thus, a scanning electron microscope (SEM) under a backscattered (BSE) mode was chosen to obtain the 2D microstructure (JEOL, Tokyo, Japan) of the aggregate with a clear outline between particles.
The simulated representative microstructures of the aggregate are used to assemble the aggregate embedded in mortar in a simulated concrete cube. XRD is then carried out, based on the thin section results that all of the minerals in the aggregate are crystalline, to obtain the mass composition of the minerals. The results are used to compare with the simulated average mineral fractions at a concrete scale.
For the porosity, simulated results are compared with MIP results. Other methods such as helium pycnometry and water absorption can be used to obtain the porosity of a porous material. However, water absorption has low liability, especially when the permeability of the material is low, such as the coarse aggregate, while Helium pycnometry results contain great errors [24].

Materials
The particle size of the siliceous limestone is between 4 mm and 20 mm. A bulk sample of about 13 kg was shipped by the producer to the lab. The different useful densities, the water absorption after 24 h, and the specific gravity were determined following the BS EN 1097-6 [25] and are given in Table 1. Test samples are subsampled from the bulk sample according to BS EN 932-1 [26] and BS EN 932-2 [27]. The specific gravity is the ratio of the surface-dried saturation density of the aggregate to water density. After determining the physical properties, the samples were stored at 105 • C until a constant weight (around 24 h) was maintained. Afterwards, the samples were kept under vacuum until the time of the MIP test. Razafinjato et al. [28,29] investigated the high-temperature behaviour of a wide range of siliceous and calcareous aggregates used in concrete. According to their results, siliceous limestone should be relatively thermally stable below 105 • C. This means that their thermal expansion does not exceed 0.25% and no micro-cracks were found on SEM images from samples heated to this temperature. Moreover, the decarbonation phase degree of calcite and transformation phase degree of quartz, which are the main compositions of the limestone, are both above 500 • C. Therefore, temperatures up to 500 • C should not have an effect on the aggregates. Thin section petrography was used to determine the mineral constituents in the limestone. A representative test sample of around 1 kg was obtained from the bulk sample conforming to BS EN 932-2 [27]. The test sample was subsequently dried at 40 • C until a constant weight was reached. Lastly, the test sample was impregnated with epoxy resin. Then, six standard petrographic thin Sections 28 × 48 mm 2 were prepared from the test sample. Thin sections were made in such a way as to guarantee the representative that each section contains at least 20 aggregates with different size distributions. Thin sections were studied in a Leica DMLP petrographic microscope (Leica, Wetzlar, Germany) in transmitted plane polarized light (PPL) and crossed polarized light (XPL) using 5×-10×-20×-50× objectives. Micrographs were acquired using the Leica digital camera (Leica, Wetzlar, Germany) at a 2560 × 1920 resolution.

Scanning Electron Microscope
In order to cover the diversity of the aggregate, five polished sections were prepared with about 40 differently sized aggregates included. The aggregate was taken from the bulk sample according to BS EN 932-2 [27]. BSE images were acquired using a JEOL scanning electron microscope (JEOL, Tokyo, Japan) equipped with an energy-dispersive X-ray (EDX) spectroscopy detector (JEOL, Tokyo, Japan) operated at 15 kV with three magnifications of 800×-1000×-2000× (with image resolution of 1280 × 960 pixels) and a Quanta FEC 450 scanning electron microscope operated at 20 kV with a magnification of 2500× (with image resolution of 2048 × 1770 pixels).

X-ray Diffraction
Three representative aggregate samples (each around 100 g) were prepared from the bulk sample according to BS EN 932-2 [27] and then were ground using a grinding machine prior to chemical analysis. Around 9 g was then taken as the test sample from each ground sample. We also added 10 wt% (1 g) Zincite (ZnO) into each test sample. Measurements were carried out by using a Thermo-Fisher-Scientific ARL X'tra diffractometer (Thermo-Fisher-Scientific, Waltham, MA, USA) equipped with a Peltier cooled detector (Thermo-Fisher-Scientific, Waltham, MA, USA), operated at 40 kV and 35 mA using CuKα, from 5 • to 75 • 2θ at 0.02 2θ increments with 1 s counting time increment. The X-ray wavelength was 1.54 A.

Mercury Intrusion Porosimetry
MIP measurements were undertaken to investigate the pore size distribution of the aggregate. In order to contain the diversity of the pore distribution, six aggregates with a size of around 4 mm were taken from the sample kept under vacuum. For each measurement, around 1.3-2.3 g of the sample was used. Intrusion and extrusion were done in a low-pressure and high-pressure PASCAL 140 Series instrument (Thermo-Fisher-Scientific, Waltham, MA, USA). First, the sample was tested in the low-pressure instrument, where the pressure was gradually increased from 0 to 200 kPa and followed by an extrusion phase, during which the pressure was decreased to 100 kPa. Then, the sample was moved to the high-pressure instrument, where the sample was subjected to mercury pressures ranging from 0.1 to 200 MPa and back to 0.1 MPa. As such, the sample went again through an intrusion and an extrusion phase. The pore size distribution data were computed based on the Washburn equation [30]. The surface tension of mercury is 0.485 N/m 2 , and the contact angle of mercury with a limestone surface is 130 • .

Framework Development of 3D Microstructure Simulation
The microstructure simulation process can be divided into three parts. Firstly, the petrographic characteristics of the limestone and its 2D microstructure from SEM-BSE images were described in Section 4.1. The image threshold and binarization process were then elaborated in Section 4.2. Finally, Section 4.3 illustrates how to simulate the 3D microstructure based on the 2D data.

Thin Section Petrography and the 2D Microstructure Features from SEM-BSE
The minerals composing the limestone were identified based on their optical properties, as illustrated in [31]. The limestone is composed of micrite (microcrystalline calcite (CaCO 3 ) with a crystal size less than 5 µm) with minor dolomite and interspersed microcrystalline and cryptocrystalline quartz (SiO 2 ), as shown in Figure 1. Quantitative results (Section 5.3.1) show that the mass contents of calcite and quartz in this limestone are around 78.7% and 13.21%, respectively. The rock also contains bioclastic fragments (fossil relics) filled or partly replaced with micro-to crypto-crystalline quartz ( Figure 1b). The micro-to crypto-crystalline quartz occurs in some domain conspicuously dispersed throughout the carbonate (Figure 1d), partly also clustered in a vein (Figure 1c). A small amount of well-crystallized sparite (calcite crystals of over 5 µm) was also discovered, as shown in Figure 1a. According to [32], the alkali reactivity of this limestone mainly attributes to microcrystalline to cryptocrystalline quartz.
Two-dimensional images with a clear outline between the particles are the prerequisite for 3D microstructure simulation. However, the Leica DMLP petrographic microscope (Leica, Wetzlar, Germany) in the laboratory is limited for this purpose, due to the pixel size limitation (200 × 200 µm 2 under a magnification of 50×). Ever since the pioneering work by Scrivener and Pratt [33], BSE imaging has been widely used as a technique to examine the microstructure of cement-based materials [34,35]. In this study, we extended it to explore the microstructure of the coarse aggregate. Two representative SEM-BSE images with three energy-dispersive X-ray spectroscopy (EDX) analyses are shown in Figure 2. In a SEM-BSE image, phases with a higher atomic number can scatter more electrons so that they are brighter, while phases of a lower atomic number absorb more electrons and are darker. The brightness of a pixel is represented by its grey value on the SEM-BSE image. The higher the grey value is, the brighter the pixel. The grey value is positively related to the backscattering coefficient η, which is the fraction of scattered electrons. Based on this principle, in order to identify different phases on the BSE image, we calculated the backscattering coefficients for the minerals that were discovered by XRD, namely calcite, dolomite, fluorite, and quartz. The calculated coefficients η for them are 0.161, 0.143, 0.186, and 0.147, respectively [36]. Based on the calculated backscattering coefficients and the EDX images, in Figure 2a, the dark-grey part is quartz, and the black part is pore, while the light-grey part is calcite. Additionally, dolomite (CaMg(CO 3 ) 2 ) (red circle in Figure 2e) shows a similar dark-grey color to quartz, as their backscattering coefficient is very close to each other (0.143 and 0.147). However, it was found that in this siliceous limestone, the dolomite shows a cuboid crystal shape (red circle in Figure 2e), while the shape of the quartz crystal is relatively irregular. By checking the crystal grains visually, most quartz grains have a size between 1 µm and 10 µm, with part of the quartz crystals clustered into larger grains. The SEM-BSE images combined with EDX analysis confirms the presence of micrite and micro-to crypto-crystalline quartz and their distribution pattern, as analyzed from the thin section results. As stated previously that the main purpose of this work was to obtain the spatial distribution of the ASR-related phases (pore and reactive silica) for further ASR simulation, we only segmented quartz and pores from the background, as will be explained in the next section. The unreactive minerals (mainly of calcite and dolomite) were viewed as one single component, and are called "the others" in the continuation of this paper.

Thresholding & Binarization
A high-quality original SEM-BSE image is the prerequisite for accurate segmentation of features and subsequent quantification steps. More than 200 SEM-BSE images were taken under the magnification of 1000×, 2000×, and 2500×. However, as described above, the similar grey value between dolomite and quartz ( Figure 2e) makes it difficult to segment quartz from dolomite. In order to avoid this issue, those SEM-BSE images that only contain quartz, pore, and calcite were selected based on the shape of the particles and the grey value, combined with EDX images. In the end, 90 SEM-BSE images were selected for binarization. Only after binarization, the pore area fraction and quartz area fraction can be calculated and their fraction distribution patterns can be analyzed. The area fraction is viewed as the volume fraction based on the 'Delesse principle' [37], and both are referred to as a fraction in the continuation of this paper. For optimum performance of binarization, images were firstly pre-processed to remove background noise through a median filter [38]. Each output pixel contains the median value in a two-by-two neighborhood around the corresponding pixel in the input image.
Image binarization means converting a grey-scale image with pixel values ranging from 0 (black) to 255 (white) in a 8-bit image into a binary image with a pixel value of either 0 (black) or 1 (white). The image will be correspondingly separated into two phases by comparing each pixel's grey value to a grey-scale threshold. Binarization methods can be distinguished by the applied method to determine the threshold. Each method has their own equation and optimal parameters, as compared and concluded by Kim [39].
Due to the low pore fraction of the limestone, a separate pore peak does not occur in the grey-scale histogram of most of images, as seen in Figure 3a. Thus, the overflow segmentation method of Wong [40] was used in order to segment pores accurately. This method calculates the inflection point (the intersection between two linear segments, as shown in Figure 3b) of the cumulative area segmented curve as the upper grey-scale threshold, where a small incremental grey value will cause a sudden increase of pore fraction. The grey-scale threshold of quartz was calculated using Ostu's standard method [41], since there is a strong contrast between quartz and the background. This method utilizes the zeroth-and the first-order cumulative moments of the grey level histogram and the discriminant analysis to calculate the threshold so that the intraclass variance of the thresholded black and white pixels is minimized. After thresholding, images were binarized for the pore and quartz, repetitively. The fractions of pore, quartz, and the others were then calculated from the binarized images. Figure 4 shows two original images with its binarization images of pore and quartz. It can be visually seen that the location, particle size, particle distribution of the pore, and the micro-to crypto-crystalline quartz have been separated from the background very successfully using the above segmentation methods.   As stated in Section 1, the REV of the coarse aggregate used in concrete is usually bigger than 1 × 1 × 1 mm 3 [15]. If we want to simulate the real 3D microstructure of the siliceous limestone in such a relatively large size, a 2D image with a size of at least 1 × 1 mm 2 would be needed. However, the microstructural information, such as the microto crypto-crystals, will not be captured in such a large image considering the pixel size limitation. Instead of trying to obtain such a large REV, we simulated a few separated representative 3D microstructures (a size of 100 × 100 × 100 µm 3 ) from 2D images at a microscale so that the real microstructure information at a microscale from 2D SEM-BSE images can be captured as much as possible. These representative 3D microstructures at a microscale (µm) are then used to assemble the real aggregate in concrete at a mesoscale (mm). For more clarity, a brief schematic diagram of this method is shown in Figure 5. Figure 5a shows a real 3D irregular coarse aggregate with a size of 4 mm. Figure 5b illustrates a random representative 2D slice of its microstructure which is composed by a few representative microstructures, as represented by different colors. An application of this method and the upscaling process is further detailed in Section 5.3.2. With this method, a REV of the aggregate is avoided. We set the simulation voxel size at a microscale as 1 × 1 × 1 µm 3 , considering the quartz particles are almost all above 1 µm. The minimum pore size was not considered because it can be as small as a few nanometers, which would cause a big burden on the computer. More random-access memory would be needed and the simulation would be very slow if the voxel size is set to be a few nanometers. In order to match the simulation voxel size, the original SEM-BSE images with a size of 1280 × 960 pixels (0.1 × 0.1 µm 2 /pixel) were then resized to a size of 128 × 96 pixels with a pixel size of 1 × 1 µm 2 and binarized again. The influence of resizing on the pore and quartz fraction and distribution is discussed in the next paragraph. Resizing was realized using the bicubic interpolation method [42] that the output pixel value is a weighted average of pixels in the nearest 4-by-4 neighborhood. One can also skip the resizing step, but must be aware that more computation resources will be needed if the pixel size is too small.
In order to investigate the influence of resizing on the mineral fraction, we calculated the pore and quartz fraction for each binarized image and the cumulative average fraction for pore and quartz before and after resizing, respectively, using all of the 90 binarized images. Figure 6 shows the results. The blue curve in Figure 6a shows that the average pore fraction of the limestone is around 4% before resizing, while the red curve shows that it reduced to around 1% after resizing. Figure 6b shows that the average quartz fraction before resizing is around 19% and becomes 21% after resizing. It can be seen that after resizing, the pore fraction was greatly reduced by more than 50%, while the quartz fraction only fluctuated in a small range (within 2%).
(a) (b) Figure 6. Cumulative average fraction of pore and quartz within 90 binarized images before and after image resizing. (a) Cumulative average pore fraction before and after image resizing; (b) Cumulative average quartz fraction before and after image resizing. The x-axis is the number of images used to calculate the cumulative average pore or quartz fraction. The y-axis is the cumulative average pore or quartz fraction. Figure 6 also shows that the minimum number of images required to give a stable average fraction is around 30 for both pore and quartz. Figure 7 shows the 2D pore and quartz distribution change after resizing of the original SEM-BSE image in Figure 4a. It can be concluded that after resizing, the pore distribution becomes much sparser, while the quartz particle distribution almost remains the same, except that a small part of the isolated particles are not withheld. This result indicates that the pore size in this limestone is at both a nanoscale and microscale, with most of the pores at a nanoscale, while the quartz crystal size is mostly at a microscale. What should be kept in mind is that the pore fraction to be simulated is actually the air void fraction (over 1 µm), since our simulation voxel size is 1 × 1 × 1 µm 3 .

Representative Image Selection
The distribution of minerals in the siliceous limestone is very heterogeneous, especially the quartz distribution. Therefore, representative 2D SEM-BSE images with the size of 128 × 96 µm 2 should be selected carefully for 3D simulation. In our work, we assumed that the mineral fraction follows a lognormal distribution in reality. It is a convenient and useful model in engineering sciences [43]. If x has a normal distribution, then the exponential function of x, y = exp(x), has a lognormal distribution. A random variable which is lognormally distributed takes only positive real values. The representative images were then selected in the following way: • Firstly, the normalized histograms of pore fraction and the quartz fraction from all of the 90 binarized images were drawn as shown in Figure 8. In Figure 8, the fraction distribution is divided into 10 bins represented by the blue bar. The x-axis is the pore or quartz fraction. The y-axis represents the probability distribution of the pore or quartz fraction within 90 binarized images. • Based on the histogram, the lognormal distributions were fitted for the pore and quartz fractions, respectively. Equation (1) is the lognormal probability density function.
where X is the pore or quartz fraction, and P is the probability density at fraction X. µ is the mean value and σ is the standard deviation. These are the mean value and standard deviation of the variable's natural logarithm, not the expectation and standard deviation of X itself. µ and σ can be calculated through µ = mean(ln(X)) and σ = std(ln(X)), respectively. The fitted (µ σ) based on the data from 90 images for pore and quartz are (1.20, 0.50) and (2.70, 0.52), respectively. The functions are drawn in Figure 8 by a black curve. Three representative statistical variables (the mode, the mean, and the median with their definition shown below Table 2) of the pore and quartz fraction calculated from the fitted lognormal distribution and images are shown in Table 2. The differences of these variables of pore fractions between the fitted curve and original data from 2D images are very small, except the mode pore fraction of the fitted curve is around 1% lower than the original one. The differences of these variables of quartz fraction between the fitted curve and raw data are also small, mainly below 2%, except that the median quartz fraction of the fitted curve is 4% higher than that of the original data. Generally speaking, the fitted lognormal distribution of the pore and quartz fraction can represent the true distributions pore and quartz fractions in the selected 90 images.  1 The mean fraction value. 2 The fraction value separating the higher half from the lower half. 3 The global maximum fraction value.

•
Finally, representative fractions were selected through the probability density function curve. In total, 10 representative images were selected: Five images with quartz fractions of (6.58%, 10.44%, 20.68%, 32.58%, 41.23%) and another five images with a pore fractions of (0.66%, 2.18%, 3.26%, 4.3%, and 8.0%) as indicated by the x-axis values of the red solid circles in Figure 8a,b, respectively. These five fractions are representative not only because they cover the whole fraction distribution span of silica or pores, but also capture the main characteristics of the fraction distribution, such as the global maximum fraction value. Figure 9 shows the correlation between the pore and quartz fraction distribution on 90 images. It can be seen that except in parts of images where the pore fraction remains relatively stable, almost half of the images (indicated by the blue part, around 43 images) have a positive correlation between the pore and quartz fractions, which means that as the quartz fraction increases or decreases, the pore fraction increases or decreases correspondingly. Therefore, the selected pore and quartz fractions are positively correlated in the simulated 3D microstructures. The final five simulated 3D microstructures (with microstructure ID from 1 to 5) have pore fractions and quartz fractions as follows: Microstructure 1: (0.11%, 6.48%); Microstructure 2: (0.51%, 10.41%); Microstructure 3: (0.93%, 20.51%); Microstructure 4: (1.45%, 32.42%); and Microstructure 5: (2.77%, 40.91%). Again, the simulated pore fraction is actually the air void fraction pore size over 1 µm.

Simulation Method
According to Quiblier [44], a 2D image of a porous medium can be characterized by two functions: (1) A probability density function (PDF) or the histogram of an image; or (2) an autocorrelation function (ACF). The physical meaning of the ACF in this paper is a measure of the probability that two pixels on an image at a specific distance belong to a same phase. Most importantly, these two characteristics remain unchanged for two and three dimensions. It should be noted that the porous medium does not need to be isotropic. An anisotropic ACF can be obtained by analysing several mutually perpendicular thin sections or SEM images.
Let I(i, j) be the value 1 or 0 depending on whether the pixel P(i, j) is located in a phase s of interest or not. The probability function for an M × N image is determined using Equation (2) [22]: where s is the target phase, which is pore or quartz in this paper. The autocorrelation function of phase s is evaluated using the following Equation (3): where x, y are the coordinate values of a pixel on the image, and 1 ≤ x ≤ M, 1 ≤ y ≤ N.
One may recognize that Equation (3) is a convolution where the kernel is rotated 180 • , and as such, can be performed rapidly using Fourier transform methods to improve the calculation. Given the two-dimensional calculated S 1 (x, y), we can obtain the desired onedimensional or isotropic autocorrelation function S 3 (r) by averaging over S 2 (r, θ) values, whose value should be calculated through the Equation (4), at a fixed radius r, as shown by Equation (5).
where the value S 1 (r cos θ, r sin θ) can be calculated by bilinear interpolation from the value of S 1 (x, y).
After getting the one-dimensional autocorrelation function S 3 (r), the 3D microstructure reconstruction can be achieved by four steps.
Step 1: Generation of a random number N(x, y, z) at each point of the target domain following a normal distribution. All N points are completely uncorrelated. The initial phase at each point in this domain is considered as the others in our work.
Step 2: For each point (x, y, z), considering a neighbouring domain composed of points (x + i, y + j, z + k) up to a certain distance which corresponds to the number of steps beyond which ACF becomes negligible (10 µm for quartz and 5 µm for pores in this paper). The values N in this neighbourhood are then linearly combined so as to form a new number denoted as R(x, y, z), which is correlated and defined by Equation (6). The periodic boundary condition is applied in the X, Y, and Z directions during the calculation of R in our work.
It has been proven by Bentz [22] that the 3D ACF F(r) calculated from the 2D ACF S 3 (r) by Equation (7) can be expressed as a function of the coefficient F(x, y, z).
where S 3 (0) is the area fraction of the target phase on the image.
Step 3: A threshold R needs to be operated on the filtered domain to match the approximate 2D fraction of each phase of interest.
Step 4: Repeating Steps 1 to 3 until all the phases of interest are reconstructed. According to Bentz [22], further modification may be required to adjust the 3D ACF through adjusting the hydraulic radius of the interested phase. However, in our method, this issue is solved in Step 2 by setting the neighbour range of each point to the 2D ACF "pessimum value", and no further adjustment is needed.

Experimental Corroboration of the Simulated 3D Microstructure
Five 100 × 100 × 100 µm 3 microstructures at a microscale with a voxel size of 1 × 1 × 1 µm 3 were simulated based on the methods described previously. Firstly, the 3D pore distribution with different fractions was simulated in each domain of 100 × 100 × 100 µm 3 based on the five 2D selected representative images of pores. After that, the 3D quartz distribution with different fractions was simulated in each domain based on the five 2D selected representative images of quartz. A periodic boundary condition in the X, Y, and Z directions was used during the simulation. Figure 10 shows the final five simulated 3D microstructures. From Microstructures 1 to 5, the pore fraction increases from 0.11% to 2.77%, and the quartz fraction increases from 6.48% to 40.91%. The pore and quartz fraction is positively correlated in each simulated 3D microstructure, as described in Section 4.3.1. Once again, the simulated pore fraction is actually the air void fraction (pore size over 1 µm).
The simulated particle sizes and their distribution characteristics of pore and quartz were first checked visually by comparing the random 2D slices selected from the 3D simulated microstructures with an original 2D image. Numerical checking was subsequently completed so that the ACF calculated from the simulated microstructures was compared with the original ACF of the original 2D images. Then, the five microstructures were upscaled to simulate the irregular coarse aggregate embedded in a 40 × 40 × 40 mm 3 concrete cube. The upscaling technique is detailed in Section 5.3.1. The average volume fraction of the quartz and air void fraction of the simulated aggregates in the concrete cube were compared with the quantitative results from XRD and MIP. The whole simulation work was implemented using MatLab-R2019b (R2019b, 2019, The Mathworks, Natick, MA, USA). Figure 10. The final five simulated 3D microstructures of the siliceous limestone from 2D images. (a) Microstructure 1 with a pore and quartz fraction of (0.11%, 6.48%); (b) Microstructure 2 with a pore and quartz fraction of (0.51%, 10.41%); (c) Microstructure 3 with a pore and quartz fraction of (0.93%, 20.51%); (d) Microstructure 4 with a pore and quartz fraction of (1.45%, 32.42%); (e) Microstructure 5 with a pore and quartz fraction of (2.77%, 40.91%). There are three phases in all microstructures: 1-pore, 2-quartz, and 3-the others. Figure 11 shows six random slices of quartz distribution from the simulated 3D microstructure 1 (lowest quartz fraction) and 5 (highest quartz fraction) and their original binarized images. Figure 11a shows that at a lower quartz fraction, the quartz particle sizes are smaller and distributed sparsely in the simulated model, which is similar to the original image. For a higher quartz fraction in Figure 11b, quartz particles are clustered together, which makes the structure more closely connected. This characteristic is also simulated relatively successfully if one compares it with the original one. Figure 12 is an example of a 2D pore distribution comparison between the 2D sections from the simulated Microstructure 5 with the highest pore fraction of 2.77% and its original image. Due to the low pore fraction above 1 µm (average is around 1%, as shown by the red curve in Figure 6a) of the limestone and the resizing effect, the pore is scattered sparsely. The pore size and distribution pattern also looks similar to the original one. Therefore, to the naked eye, a similarity undoubtedly exists in the size and 'scattering' of the pore and quartz between the simulated microstructure and the real image.

Numerical Checking
The ACF curves of pore and quartz in the original 2D images are the input for the 3D microstructure simulation. In order to check if the simulation successfully remains as the characteristics of the original ACF, in this section we computed the ACF curves from the five simulated 3D microstructures and compared them to those of the input original images. Figure 13 shows the ACF curve of quartz. The red dashed curve is the ACF from the simulated 3D microstructure, while the black curve is the ACF from the original 2D image. It is obvious that the original and reconstructed ACF are seen to nearly overlap in the distance range of a value from 0 to around 10 µm for quartz in all of the subfigures. As stated before, the ACF describes the probability of two random points with a distance r located in the same phase. The distance range of 0 to 10 µm is actually the PSD of quartz, and 10 µm is the maximum 2D particle size of quartz. Therefore, we can conclude that the simulated PSD agrees very well with the original PSD. At a longer distance, the curves do not overlap too well. The original ACF curve fluctuates due to the random nature, as also found in [22], while the simulated ACF curve is flat due to the extent of the filter operation used in the generation algorithm. However, this small difference can be ignored, since it does not influence the PSD.
Except the PSD, the ACF curve contains more information, such as the mineral fraction α A , the value φ A 2 when the ACF function becomes stable, the characteristic particle size L cd , and the specific surface area α sur f (ratio of the phase surface area to the total volume of the sample) of the phase of interest [20,45].
Taking Figure 13a as an example, the value of the curve at r = 0 corresponds to the 2D quartz fraction from the 2D image (black curve) and the 3D quartz fraction from the corresponding reconstructed microstructure (red dashed curve). The plateaus value of the ACF function corresponds to φ A 2 . The intersection of the slope of ACF at r = 0 and the plateau region indicates a characteristic dimension of the quartz structure l cd . This parameter can be used to calculate the average quartz particle size d av according to [46]. The slope of the ACF at r = 0 relates to the specific surface area α sur f of the quartz. By comparing the ACF curves in Figure 13, we can say that all the mentioned parameters of the simulated microstructures match very well with that of the original images, except Figure 13c, where the plateaus value of the original image is a little bit lower than the simulated one.    Figure 14 shows the corresponding ACF curves of the pore. The simulated ACF and original ACF almost overlap in the distance range from 0 to 5 µm, except in Figure 14c with a pore fraction of 0.93%, where the simulated PSD is narrower than the original figure. However, the pore fraction, the slope at r = 0 of the red dashed and black curves in each subfigure are the same, as well as the plateaus value. Generally, the same conclusion as for quartz can be made for pores, that the simulated ACF of pores agrees well with the original one.

Quantitative Results from XRD
The XRD pattern in Figure 15 confirms the observations from the thin section and SEM results that this limestone is composed mainly of calcite (volume fraction of 78.7%) and quartz (volume fraction of 13.21%) with a small part of dolomite and fluorite (volume fraction of 8.08%). The mass and volume fractions of the minerals are shown in Table 3. The volume fractions were calculated by dividing the mass fraction (obtained using the software TOPAS (Academic V4.1, 2007, Coelho Software, Brisbane, Australia) by the corresponding density. The density of calcite, dolomite, fluorite, and quartz is 2.71 g·cm −3 , 2.84 g·cm −3 , 3.13 g·cm −3 and 2.65 g·cm −3 , respectively, as obtained from [47].  2 1 The average volume fraction of all phases other than quartz. 2 The average volume fraction of the quartz.

Upscale from Microscale to Mesoscale
In order to compare the simulated mineral (pore, quartz, the others) volume fraction with the real volume fraction from XRD or MIP, a concrete cube of 40 × 40 × 40 mm 3 was simulated. The size of 40 × 40 × 40 mm 3 was chosen because the minimum size of a REV of concrete at a mesoscale should be at least 2.5 times bigger than the maximum aggregate size (16 mm in this paper) [48,49]. The cube is composed of mortar and irregular aggregate with a particle size from 4 mm to 16 mm. An aggregate volume fraction in concrete of 30% was used for simulation. The microstructure of the coarse aggregate was composed by the five representative simulated 3D microstructures. The simulation process and upscaling technique is described below.
Firstly, the voxel size of the concrete cube was set to be 100 × 100 × 100 µm 3 , which means that the cube size of the simulated 3D microstructure is now the voxel size at mesoscale. Therefore, the final simulated cube has a lattice size of 400 × 400 × 400.
Then, the irregular aggregates from 4 mm to 16 mm were distributed randomly inside the cube using the ANM model of Qian [50]. The volume fraction of the aggregate at different sizes from 4 mm to 16 mm was calculated based on the cumulative distribution function (CDF) of the coarse aggregate, as shown in Figure 16. The surface of the simulated cube and a random slice of the simulated cube are shown in Figure 17a,b. After this step, each node in the cube is 0 or 1, where 0 represents mortar and 1 represents the aggregate. By accounting the node number of 1 of the cube, we found that the total nodes occupied by aggregates were 19,136,000 and the simulated aggregate volume fraction in the cube was 29.9%, which is very close to the true volume fraction.
The last step was to allocate the five simulated 3D microstructures to these 19,136,000 nodes. Numbers 1 to 5 were used to represent the five simulated 3D microstructures, respectively, as shown in the first column in Table 4. The volume fraction occupied by each microstructure was calculated according to the fitted lognormal distribution function of the quartz fraction, as determined in Section 4.3.1 (Figure 8b, X∼Lognormal (2.70, 0.52 2 )). The calculated volume fraction of each microstructure was shown in the third column in Table 4. Then, the number of nodes for each microstructure was calculated to be 1.053 × 10 6 , 1.2967 × 10 7 , 3.862 × 10 6 , 7.90 × 10 5 and 4.96 × 10 5 , respectively, as shown in the second column in Table 4. Subsequently, we randomly assigned 1.2967 × 10 7 nodes out of node 1 to be 2 (microstructure 2) and repeated this process until each node 1 was reassigned. A schematic diagram of the final microstructure of one aggregate at mesoscale is shown in Figure 17c.   Table 4 also shows the simulated average mineral fractions (mineral volume/aggregate volume) of the concrete cube. The simulated average quartz fraction at mesoscale is around 13.93% (Table 4, Column 5). It only increased a factor of about 1.05 compared with the real quartz fraction (13.21%) obtained from XRD. Consequently, the fraction of the others reduced from 86.79% to 85.43%. It can be seen that the simulated compositions agrees very well with the XRD results only with a relative error less that 6% for both the quartz and the others. Figure 18 shows two representative pore size distribution patterns of the limestone: a low pore fraction of 0.23% (Figure 18a) and a high pore fraction of 2.56% (Figure 18b) obtained from MIP. It can be seen from the blue curve in Figure 18 that the pore size from the MIP results is mainly located between 10 nm and 100 nm, and over 10 µm. Most of the pores are over 10 µm, while almost no pore between 100 nm and 10 µm was found. One should remember that in our model, the simulated pore fraction only considers the aid void (pore size over 1 µm). Therefore, a pore fraction over 1 µm was extracted from the MIP results and concluded in Table 5. The results were then compared with our simulation result in Table 4. The data in the second column of Table 4 shows that the air void fraction in the simulated concrete cube mainly varies from 0.11% to 2.77%, which is close to the distribution range of 0.09% to 2.35% from the MIP results, as shown in Table 5 (column 3). Additionally, the simulated average air void fraction of the concrete cube (Table 4, column 4) is 0.65%, which also agrees well with that of MIP results (0.7%). Moreover, the average air void fraction obtained from SEM-BSE images was around 1% (the stable cumulative average value of the red curve in Figure 6a), where only 0.35% was lost during the simulation. The above analysis indicates that the simulated pore fraction distribution matches well with the results derived with MIP, and the simulated average pore fraction over 1 µm is also close to that of the MIP results, as well as of the SEM results. However, one should be aware that MIP may underestimate the true pore fraction because the mercury is unable to access all of the inner open pores or if less open pores' frequencies are on the surface, while the image analysis could overestimate the true pore fraction since it is the 2D intersection of the 3D pore being analyzed, and the binarization process can induce error as well.

Application of the Model
The main purpose of obtaining the microstructure of the alkali-reactive aggregate is to use it as a basis in a ASR multiscale simulation model. Such a model starts the simulation of ASR from the chemical reaction at a low scale to the physical response both at low and high scales. Each of the five separated simulated microstructures of the siliceous limestone can be used as the domain of interest at a microscale to enable deep exploration of the chemical reaction mechanism of ASR based on the methodologies we proposed in [51]. For example, the chemical reaction sequence, the location of the ASR products or the initial expansion sites (inside or outside the aggregate), or the influence of silica fraction on ASR, and so forth, can be clarified as realistically as possible. Moreover, the fracture process at microscale (µm) can be simulated once the information about the ASR products are obtained. Based on the multiscale fracture modelling technique in [50], by passing the mechanical simulation results at microscale to mesoscale, the cracking process of the concrete cube can be simulated then. In this way, the chemical reaction and cracking mechanism induced by ASR at different scales can be captured as much as possible so that better suggestions for preventing ASRs in newly-built structures can be put forward.

Conclusions
This paper elaborated how to simulate a realistic 3D microstructure of a siliceous limestone based on the autocorrelation function obtained from image analysis performed on 2D SEM-BSE images of the aggregate. The following concluding marks can be made.
(1). Suitable experimental methods should be chosen to obtain 2D images of the aggregate with a clear outline between particles, considering the particle sizes of the target phase, and suitable binarization methods should be chosen to segment the target phase from the images. (2). The pore and silica fraction distributions on the 2D SEM-BSE images of the limestone can be fitted by a lognormal distribution, which can be used to select representative images for 3D microstructure simulations. (3). The simulated microstructures of the siliceous limestone are able to retain the visual characteristics, such as the particle shape and spatial scattering, as well as the statistical characteristics of the air void and quartz silica, such as the fraction, characteristic particle size, and the specific surface area of the parent limestone. However, the pore information below the simulation voxel size (1 µm) is lost during the simulation. (4). The simulated microstructures can be used to assemble the aggregate at a mesoscale embedded in mortar following the obtained lognormal distribution. The average air void and silica fraction show good consistency with the XRD and MIP results. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data is contained within the article.