Simulation of Permeation of Saturated Cement Paste Based on a New Meso-scale Pore Network Model

This paper presents the simulation of the permeation of saturated cement paste based on a novel pore network model. First, a 2D hydration model of cement particles was developed by extending the work of Zheng et al. 2005 to provide the background for the network construction. Secondly, the establishment of the pore network model and simulation of permeation of saturated cement paste were carried out. The irregular pores between any two hydrated cement particles were linearized with clear distances as the diameters of pores. The straight tubular pores were interconnected with one another to form the network model. During this process, the weighted Voronoi diagram was employed to operate on the graphical expression of the hydrated cement particles. Water permeation in saturated cement paste was simulated to verify the pore network model. Finally, the factors including water–cement ratio, reaction temperature, reaction time and cement particle size that would influence water permeation were numerically investigated.


Introduction
Cementitious materials have been extensively used in many fields of construction because of their low price, abundant raw materials, high compressive strength and relatively simple production and construction techniques. However, it is difficult to ensure that the cementitious materials exposed to natural conditions meet their design durability requirements. The durability problem is related to specific material intrusion such as water, chloride ions, sulfate ions, oxygen, and carbon dioxide. Therefore, the study of permeability of cementitious materials provides an important approach for solving the durability problem.
Page 2 of 15 Zhou et al. Int J Concr Struct Mater (2021) 15:36 properties of cement-based materials including autogenous shrinkage of cement pastes (Koenders, 1997) and transport properties of cement pastes (Ye et al., 2006). CEMHYD3D is another 3D microstructure-based simulation model for cement hydration developed by Garboczi (Bentz, 1997, 2005;Bentz & Garboczi, 1991;Garboczi & Bentz, 1992). CEMHYD3D is a latticebased approach based on digital image technology and considers the effect of the reaction stoichiometry of the cementitious material. Numerical modeling approaches are widely used to establish the link between the cement hydration microstructure and its permeability. Among the different numerical methods, Lattice Boltzmann method is one of the most utilized method to predict the permeability of cement paste. Zalzale and McDonald (2012) have used Lattice Boltzmann method to investigate the permeability of cement paste using the numerical models CEMHYD3D and µIC (Bishnoi & Scrivener, 2009). It was concluded that the simulated permeability is larger than the experimental data, while the results obtained from µIC compare favorably with CEMHYD3D. Zhang (2017) presented a Lattice Boltzmann modeling approach for estimating the relative permeability of cementitious material, where the 3D microstructure was obtained from high-resolution X-ray micro-CT. The good agreement of the simulated result with experimental result indicates that the moisture distribution and permeability of cement paste are strongly dependent on its microstructure. Besides, other numerical methods have been developed to estimate the permeability of cement paste. Sun et al. (2014) numerically calculated the transport properties of a microscale cement paste specimen by applying the transmission X-ray microscope (TXM) characterization techniques and a finite difference method based computational program. Yu et al. (2018) presented a fractal model to predict the transport properties of cement paste. Furthermore, Li and Xu (2019) proposed a microstructural-based numerical test method based on three-dimensional (3D) finite element method (FEM) to describe the pore-scale flow and compute the permeability of hydrating cement paste, where HYMO-STRUC3D model was used for the microstructure evolution of cement paste.
Another frequently used tools for the numerical simulation of permeability in porous media is the network model. In 1956, Fatt (1956 proposed a network model using resistors as analogs of porous medium which could effectively reflect the microstructure of porous materials. The network models when combined with a simple percolation theory can quantitatively describe the macroscopic properties such as the permeability, relative permeability, capillary pressure of materials (Larson et al., 1981) (Blunt et al., 1992). The advantage of the network model is that it can explain macroscopic behavior by explicitly analyzing the applicable physics at the pore level. However, when using the network model, the internal structure of the material must be accurately analyzed and constructed to make the calculation result more accurate (Arns, 2004). Ye et al. (2003) used HYMOSTRUC3D model to develop the hydration model and subsequently, converted it into a linked list pore network model to determine the unidirectional absolute permeability of cement paste (Ye et al., 2006). Stroeven et al. (2012) have developed a novel approach to the pore network in virtual concrete using Random Node Structuring algorithm. The algorithm involves the generation of nodes uniformly at random in model space followed by the elimination of nodes located in solid phase while connecting the remaining nodes to form cluster of capillary pores for the pore structure analysis.
A simple 2D pore network model was proposed in this paper. It started with the formation of particle. According to the size distribution of cement particles combined with the Monte Carlo method, cement particles were generated and hydrated based on the hydration kinetics referring to the study of Zheng et al. (2005b) and Wu (2015). Then the hydrated model was converted into the pore network model using the weighted Voronoi diagram to analyze the permeation of water in cement paste. There are two reasons to use the weighted Voronoi diagram. First, when the modeling space is tessellated, the size of the particles is considered, which ensures the interface, i.e., the transport path between any two adjacent particles will not intersect with any particle. Compared with the network based on the Voronoi diagram, which simply bisects the line linking two adjacent nodes, the model based on the weighted Voronoi diagram is more physical and more reasonable. Second, the interfaces produced by the weighted Voronoi diagram are planes (3D) or straight lines (2D). This keeps the advantage of simplicity of the Voronoi diagram-based network model. To the best of the authors' knowledge, the weighted Voronoi diagram has not been used in past to construct the applicable network model based on randomly positioned hydrated cement particles. The good match of the numerical results and the existing experimental data proves the validity of the proposed network model. The computing program used in this study is written in MATLAB. The proposed model is much simpler and faster than the existing hydration models like HYMOSTRUC3D, CEM-HYD3D, µIC and can also give the acceptable results. These existing models can predict the hydrated product, which is not necessary when the concern is just the water transport problem.

Model Size Selection
The model size should be consistent with the representative volume element (RVE) size of cementitious materials. Thus, it should accommodate enough cement particles to reflect the physical characteristics of hydration reaction, but should be minimized to increase the model efficiency. The model size in 2D space with each side length L , is initially set to 200 µm × 200 µm and will be discussed in Sect. 4.4. The simulation was carried out for 100 times where each 2D simulation is like a CT slice and it has been assumed that the average result of those 100 simulations approximate the 3D behavior of the real pores.
In HYMOSTRUC model, the cement particles are modeled as randomly distributed digitized spheres and the hydrating cement grains are simulated as growing spheres (Koenders, 1997;Van Breugel, 1991). However, the influence of particle shape on the cement hydration process should be considered. It has been observed that the degree of hydration is more pronounced at early ages for real-shape particles than for spherical particles (Bullard & Garboczi, 2006). Since the focus of this paper is the transport properties rather than the early ages of hydration process and considering the computational time and efficiency of the model, the cement particles are modeled as spherical in shape. The cement particle size distribution for Cement 115 and 116 issued by ASTM-sponsored Cement and Concrete Reference Laboratory ranges from 3 to 37 µm (Bentz, 1997). In the present simulation, minimum diameter ( D min ) and maximum diameter ( D max ) of the particles are taken as 2 and 40 µm, respectively. Initially, the model is assumed to be filled with two phases of water and cement particles. Thus, we have From above, the total volume of cement particles V c is: where m w , m c , ρ w and ρ c are the total mass and density of water and cement, respectively, and w/c is the water-cement ratio. The density of cement ranges 3000-3200 kg/m 3 (Wang, 2013) and taken as 3150 kg/m 3 ( van Breugel, 1992van Breugel, , 1995 in this simulation. (1)

Cement Particle Formation
The cumulative distribution function of cement particles P(D) , with diameter D are usually described by the Rosin-Rammler distribution (van Breugel, 1992) and expressed as: where α and β are empirical parameters, usually taken as 0.038 and 0.98, respectively (Zheng et al., 2005a(Zheng et al., , 2005b. Based on the stereological principles, Zheng et al. (2005aZheng et al. ( , 2005b have derived 2D cumulative distribution function of cement particles from Rosin-Rammler distribution and given as: If the probability of the ith particle P Ni , is taken as a random number within [0, 1], Eq. (4) gives the diameter of the corresponding cement particle. In the process of particle formation, the total volume of cement particles accumulated is: If V i−1 <V c and V i >V c , the total volume of cement particles generated exceeds the total volume for a given w/c ratio. Thus, the diameter of the last cement particle is: Given the model size and water-cement ratio, the number and particle sizes are generated inside the model.

Cement Particle Placement
The cement particles are randomly placed into the model by the Monte Carlo method with the condition that the center-to-center distance between any two particles is not less than the sum of their radii. The cement particle should be sorted in the descending order of the diameter before being placed into the model to facilitate the placement.

Cement Hydration Simulation
The hydration reaction of Portland cement is very complicated (Bentz, 1997). This paper follows the basic assumption for the cement hydration reaction mentioned in the literature (Stroeven, 1999).

Hydration Rate of Cement
Fig . 1 shows the hydration profile of a cement particle, where r in , r out and r air represent the radius of the unhydrated cement particle core, hydrated gel surface and air layer, respectively (Stroeven, 1999). The thickness of hydrated gel layer is given by: The hydration process of cement is believed to be dominated by two mechanisms. The first mechanism is controlled by the interface between the cement particles and water ( Van Breugel, 1991) and its reaction rate K 0 is expressed as: When the thickness of hydration product reaches the critical value δ tr , the reaction gradually transforms into the second mechanism, controlled by diffusion. Its corresponding reaction rate is: where b is the gel diffusion adjustment constant, K T is the reaction rate of the cement particles at the absolute temperature T , and expressed as: where E is the activation energy of the reaction; R is the gas universal constant, E/R is about 5364 K −1 ; K 293 is the reaction rate of the cement at an absolute temperature of 293 K. In this paper, the average value of K 293 is 0.1512 µm /h and b is 1.9614.

Radius of the Hydrated Gel Surface
To ensure that the simulation is close to the real situation, the change in particle size and water consumption are calculated based on the simulation in three-dimensional space. The relationship between the volume of hydrated cement, v c , and the volume of hydration product, v g , can be given by: where k 0 is taken as 2.2 ( van Breugel, 1992).
Equations (9) and (10) give the change in the radius of cement particle ( r in ) over a certain time. The corresponding volume change, v c , that completes the hydration reaction is: The radius of the hydrated gel surface can be expressed as: The clear distance between any two adjacent hydrating particles with their center-to-center distance d 0 and outer edge radii r 1,out and r 2,out , respectively, is:

Particle Interference and Degree of Hydration
The physical interference between adjacent hydrating cement particles as shown in Fig. 2 affects their further hydration, thus, the model assumes that the contact surface stops the reaction, whereas the non-contact surface continues.
The non-contact solid volume V i,s of the i th particle at any time is calculated by subtracting the hatched area from the total volume of the i th particle: where α j is the half-interference influence angle on particle i that can be calculated by the cosine theorem: The total degree of hydration of the model α c (t) at any time t can be expressed as: Page 5 of 15 Zhou et al. Int J Concr Struct Mater (2021) 15:36 where V cc is the total spherical volume of the cement particles. The maximum degree of hydration when the water cement ratio is relatively low is given by (Beton-Kalender, 2010): The total hydration reaction time in this paper is 28 days with 1 h as the basic unit. The algorithm of the hydration reaction is shown in Fig. 3. Fig. 4 shows the comparison of the present hydration simulation result with the hydration test data of Danielsson (1960), degree of hydration calculated by Ye (2003) measuring the heat released from isothermal calorimeter measurements, numerical simulation result of van Breugel (van Breugel, 1995) under the conditions of ambient temperature of 20 °C and water-cement ratio of 0.3 (Fig. 4a) and 0.4 (Fig. 4b), respectively. It can be observed that the results of different test are close, and the fitting degree is high. In Fig. 4a, the simulation result of this paper is slightly higher than the test results. This is possibly because some water is trapped among hydration products and cannot be used by other cement particles for the experiments under the condition of low waterto-cement ratio. The hydration model has been validated which provides the platform for establishing pore network model. Since the paper is not focused on the hydration, the detailed study on its microstructure, pore structure and porosity is out of scope of this paper.

Network Construction and Permeation Simulation
The network model is an ideal method for describing the spatial state and topology of complex pores. The network model can effectively reflect the microstructural characteristics of porous materials and can be quantitatively described by a simple permeation theory (Arns et al., 2004). However, it is still a challenge to determine the geometry and parameters of the network model used to study infiltration problem (Ye et al., 2006). Most of the existing methods are adapted to measure the pore parameters by adjusting the model parameters (Marchand & Gerard, 1997), but these methods are not able to accurately reflect the true microstructure of the cementbased materials. Therefore, the ideal network model of cement-based materials should be directly based on the actual or simulated microstructure formed by the cement hydration process and as simple as possible. Based on this concept, the hydration simulation model developed in Sect. 2 has been transformed into the network model using the principle of weighted Voronoi diagram. The core idea is to find the nodes in the hydration model and connect all the adjacent nodes to construct a network model. The detailed technical route, the algorithm of network model generation as well as the permeation simulation have been discussed in this section.

Weighted Voronoi Diagram
Weighted Voronoi diagram is a special case of Voronoi diagram in which the function to describe a Voronoi cell is defined by a weighted (multiplicative) distance function rather than an equidistance function. Referring to the mathematical definition of the weighted Voronoi diagram (Dong, 2009), it is proposed that the node among any three related particles A i , A j and A k should satisfy: where Rf m is the influence range of the circle O m with center p m (a possible node position) as shown in Fig. 5. d(p m , A i ) , d(p m , A j ) and d(p m , A k ) are the center-tocenter distances between circle O m and particle i, j and k , respectively. R i , R j and R k are the radii of particle i, j and k , respectively. When a new particle A n with a radius of R n is added to the group of particles, the center of particle A n is first connected to the possible node p m as shown in Fig. 5. The length of this line segment is divided by R n . p m is a node if the particle A n is not influenced by circle O m , i.e.: Otherwise, p m is not a node and needs to be deleted. The node position of the area is recalculated until all the particles in the area satisfy Eqs. (20) and (21).
The search and calculation method of this node is close to the indirect method of Voronoi diagram vector generation method. Therefore, it is necessary to establish a generalized Delaunay triangulation first.

Construction of Generalized Delaunay Triangulation
The construction of the generalized Delaunay triangulation generally adopts the point-by-point placement method. The basic idea of this method was proposed by many scholars: Lewis & Robinson (1978), Lee & Schachter(1980), Bowyer (1981), Lee & Lin (1986). The steps for constructing the generalized Delaunay triangulation are shown in Fig. 6.

Network Node Connection
In the weighted Voronoi diagram each site element has its own weight thus, the convex polygon boundaries (the transport paths) cannot be determined simply by the vertical bisector. In fact, these boundaries are obtained by connecting the network node of each triangle to the network node of all triangles with which it has a common edge as shown in Fig. 7. Finally, the network node connection is completed.

Model Boundary Processing
The cyclic boundary condition must be adopted to prevent the loss of the part of hydrated cement particle that would extend outside the model boundary.
The hydration model should be copied and distributed on each side of the model to develop the extended hydration model, as shown in Fig. 8. After the replication and translation work, the generalized Delaunay triangulation and node connection algorithm explained in Sects. 3.2 and 3.3, respectively, are applied on the extended hydration model to obtain the weighted Voronoi diagram. Once the weighted Voronoi diagram is obtained, the redundant nodes and paths outside the scope of original hydration model need to be identified and cleared.

Permeation Simulation of Pore Network
For the permeability calculation of a single pore tube, the length and diameter of the tube/transport path are first determined. The length l ij is: where ( x i ,y i ) and ( x j , y j ) are the coordinates of the first and last nodes of the tube, respectively.
The diffusion and precipitation of the hydrated products may influence the pipe diameter, so the cross-section of all the elements connecting to one node should be determined according to Yip et al. (2005). But for the simplification, it is neglected in this paper. And the pipe   diameter has been taken as the minimum of the pore width at the ends of the tube (Ye et al., 2006) and calculated according to Eqs. (14) and (15). However, the pipe diameter d ij needs to be adjusted since the pore path and the center-to-center distance of the particles on both sides are not perpendicular: where d is the distance between particles, and a and b are the direction vectors of the path and the line between the particles, respectively. As proposed by Dullien (2012), the water conductivity coefficient g ij between any two nodes i and j in the model can be expressed as: It is assumed that the network model is saturated and filled with a single liquid (water). The liquid flows in from a single boundary and flows out from its opposite boundary, and the remaining boundaries are completely impermeable. For the saturated pore network system under the action of pressure difference p ij between nodes, Hagen-Poiseuille law (Liu, 2016) gives: where Q ij is the water flow velocity between nodes i and j , the unit is m 3 /s; µ is the fluid viscosity. The fluid viscosity of water is related to the ambient temperature (Al-Shemmeri, 2012) and expressed as: where T is the ambient temperature, the unit is K, B is taken as 247.8 K; C is taken as 140 K. When the ambient temperature is 20° C, the fluid viscosity of water is 1.002 × 10 -3 Pa·s.
Using Eq. (25), the total traffic Q i at node i can be written as: For incompressible fluid, the total volume of liquid flowing in and out of a node inside the model per unit time is equal, thus, the corresponding node traffic Q i = 0, but the node pressure p i is unknown. For the nodes on inflow or outflow boundary surface, Q i is unknown with known input or output pressure p i . For all nodes, Eq. (27) can be written as: where P is node pressure vector and Q is node flow vector. The coefficient matrix G 0 is: During the permeation process, fluid flows from one node to another node, thus, g ii = 0. Equation (28) is a system of linear equations with n unknowns distributed in vector P and vector Q . Since there are many 0 elements in the coefficient matrix G 0 , the solution of the equations should be combined with the shifting partition and the pseudo-inverse matrix. Then, the total water flow Q 0 of the network is calculated by adding the Q i of the nodes that lie in inflow or outflow surfaces and given as: Finally, the total permeation rate of the model is calculated using the Darcy's law: where L is the length of the fluid in the flow direction, A is the cross-sectional area through which the fluid flows.

Validation
During the simulation, the network of the model might be in connected, completely closed, partially closed, or blockage state. To ensure the simulation is close to real situation and to minimize the interference caused by randomness of the model, permeation simulation for each set of parameters is calculated for 100 times and averaged. The model connectivity η for each set of parameters is defined as follows: where n l is the number of models in the connected state and n A is the total number of repeated calculations. Phung et al. (2013) and Jin (2015) have carried out experiment on permeability of water in cement paste using stable seepage method. The present simulation was also carried out at the temperature of 20 °C and curing age

Permeability vs. Water-Cement Ratio
Page 10 of 15 Zhou et al. Int J Concr Struct Mater (2021) 15:36 of 28 days. The permeability was determined for different water-cement ratios and the final averaged simulation results shown in Table 1 are compared with the experimental results. The simulation result is observed to be consistent with the experimental result as shown in Fig. 9. Thus, the model has high degree of confidence. The slight differences in the simulation results and experimental data can be attributed to some extrinsic factors related to experimental conditions such as the applied pressure, size of sample and the fineness of the cement used.
The influence of water-cement ratio on permeability of cement paste can be explained as: from Eq. (3), the total volume of cement particles (V c ) in cement slurry decreases rapidly when the water-cement ratio (w/c) increases. This decreases the number of cement particles in the model and increases the average distance between them. The total pore volume and connectivity between the pores will increase (also observed in Table 1). As a result, the average hydraulic conductivity of each path as well as the overall permeability coefficient of the model increases.

Permeability vs. Hydration Reaction Time
Equations (9) and (10) support that the hydration reaction time of cement particles directly determines the end point of the model hydration reaction. As the hydration reaction progresses, hydration product replaces the space originally occupied by an cement particle and accumulates on its outer edge. This expansion of the particle size reduces the average clear distance between them, thus decreases the pore connectivity as well as the permeability coefficient. Table 2 also shows that the model connectivity, average permeability coefficient of the cement paste and its standard deviation decreases with the increase in hydration time.
The simulation result is validated with the previous experimental results carried out by Ye (2003) and Banthia & Mindess(1989) and is shown in Fig. 10. Ye (2003) determined the water permeability of the cement paste using direct test method under the curing temperature of 20 °C and w/c = 0.5 and Banthia & Mindess(1989) determined the permeability of cement paste (Cement Type-ASTM III, Specimen type M3) with w/c ratio = 0.5 under normal laboratory conditions of temperature and pressure. The simulated permeability of the model is 1.5-2 order of magnitude lower than the experimental results in the  Page 11 of 15 Zhou et al. Int J Concr Struct Mater (2021) 15:36 earlier stage and is closer to the experimental results in the later stage (i.e., 28 days).

Investigation of Factors Affecting the Overall Permeability of the Model
The water-cement ratio (explained in Sect. 4), hydration reaction time (explained in Sect. 4), hydration reaction temperature (T), maximum size of cement particle ( D max ) and minimum size of cement particle ( D min ) are the major factors affecting the overall permeability of the cement paste. The control variable method has been used to analyze these factors and their influence mechanism is shown in Fig. 11.

Hydration Reaction Temperature
Increasing the early curing temperature for concrete enabled to achieve the lower long-term permeability for 6 months and beyond at 28 days (Ozyildirim, 1998). According to Eq. (11), the hydration reaction rate increases with the increase in ambient temperature.
Higher reaction temperature increases the degree of hydration of cement particles as well as the volume of hydrated product and thus decreases the permeability of the cement paste. As observed in Table 3, the permeability coefficient of cement paste decreases rapidly in the range of 10-20 °C and decreases slowly in the range of 20-40 °C. Thus, as compared with other influencing factors, the change of hydration reaction temperature has the most significant effect on the permeability of cement paste. The fitting result is shown Fig. 12.

Minimum Size of the Particle
The minimum size of the particle is the main controlling factor for a small group of cement particles that fills the gap between large particles in the model. Increasing the minimum size of particle weakens its filling capacity to the pores. This gradually increases the voids and consequently the permeability of the cement paste. Table 4 also shows that the model connectivity, the permeability coefficient, and its standard deviation increases slowly with the increase of the minimum size of cement particle. Similar influence of the particle size distribution on the permeability of cement paste can be observed in the study of Pignat et al. (Pignat et al., 2005). The fitting result is shown in Fig. 13. Fig. 11 Influence mechanism of the model permeability.

Maximum Size of the Particle
The maximum size of the particle is the main controlling factor for the large particles in cement paste which play a significant role in constructing the basic structure of the permeation network. For a given model size and water-cement ratio, if the diameter of maximum cement particle is larger, lesser particles are placed in the model which increases the voids as well as the permeability of the model. Therefore, changing the maximum size of particle will directly affect the permeability of cement paste. Table 5 also shows the increasing trend of pore connectivity, permeability coefficient and its standard deviation with increase in maximum size of particle. The fitting result is shown in Fig. 14.
It should be noted that the large coefficient of variability (> 100%) observed through Tables 1, 2, 3, 4, and 5 is mainly due to the inherent random porous medium.

Model Size
In addition to hydration reaction time, hydration reaction temperature, water-cement ratio and cement particle size, the size of the hydration model will also affect the overall permeability coefficient and the model stability. The selection of model size is critical in the hydration and permeation simulation since the model size should be consistent with the RVE size of cementitious materials.   Fig. 14 Fitting curve of upper limit of cement particle size and water permeability coefficient.
Page 13 of 15 Zhou et al. Int J Concr Struct Mater (2021) 15:36 For the quantitative analysis, the simulation was performed on model of size 100 µm , 125 µm , 150 µm , 175 µm , 200 µm , 225 µm , and 250 µm with water-cement ratio of 0.5, the maximum and minimum size of the cement particle as 40 µm and 2 µm , respectively, and the hydration reaction time of 28 days at the temperature of 20 °C. The influence of model size on its overall permeability coefficient and the standard deviation are plotted in Figs. 15 and 16, respectively.
As observed in Fig. 15, the model permeability coefficient decreases significantly with the increase in model size when the size is less than 200 µm , but when the model size is larger than 200 µm , the permeability coefficient tends to be stable. Furthermore, as observed in Fig. 16, when the model size is less than 175 µm , the standard deviation of model permeability coefficient decreases significantly with the increase of model size, but it tends to be stable when the model size is larger than 175 µm . The reason for this can be explained as: when the model size is small, the total number of cement particles in the model is small. They are likely to gather in certain area of the model which increases its overall permeability coefficient. With the increase of the model size, the number of particles in the model increases continuously and the internal permeation network becomes complicated which reduces the concentration of particles in certain area. This weakens the influence of large-size pores in the model, thus reduces the permeability coefficient. When the model size is large enough, even if some cement particles are concentrated in certain area, its influence on the overall permeability coefficient of the model is small or negligible. Thus, the overall permeability coefficient of the model will no longer change with the size of the model.
Under the current conditions, considering the stability of the model results, the model size should not be less than 175 µm and considering its accuracy, the model size should not be less than 200 µm . Thus, the optimal size of the model in this paper is adopted as 200 µm and the recommended size of the model in simulation process for the cement paste should be 5 times the maximum size of the cement particles.

Conclusion
A new meso-scale pore network model for the permeation simulation of saturated cement paste using the weighted Voronoi diagram has been developed. The weighted Voronoi diagram is used because it considers the different particle size and can produce a straight interface in the model. The minimum gap between two hydrated cement particles has been considered as the diameter of the capillary tube and this method is much simpler with enough accuracy compared with the other pore network model, which generally follows the exact profile of the pores. However, the accurate curvilinear profile of the pores is still an assumption and is only critical in hysteresis simulation, such as drying-wetting cycles, which is beyond the scope of this paper. Referring the study of Zheng et al. (2005aZheng et al. ( , 2005b, hydration simulation was completed which provides the background for establishing the pore network. The hydration simulation and the predicted permeability results show good match with the previous experimental results. The parametric study for the factors influencing permeation simulation including water-cement ratio, hydration reaction age, hydration reaction temperature and cement particle  Page 14 of 15 Zhou et al. Int J Concr Struct Mater (2021) 15:36 size has been studied. For the cement particle size, upper limit of the particle has higher influence as compared to lower limit. From the perspectives of accuracy, stability and calculation duration of simulation, the size of the model is recommended to be five times the size of maximum cement particle.