## 3 Finite Element Modeling  
The multiphase 2D benchmark model in this study is discretized using the Abaqus/CAE Mesh module. The dimensions, applied boundary conditions, and discretization strategies are shown in Figure 4. An important aspect of the meshing process is ensuring continuity between the element surfaces, which must be maintained throughout the entire meshing process in CZM simulations. A four-node quadrilateral element (CPS4) is used for both the CZM and PFM simulations by implementing a free meshing technique. For all the simulations, the 2D benchmark model is discretized with an average element length of 0.05 mm. For CZM simulations, once the model is meshed, in-house developed Python algorithms process the FE mesh file to generate zero-thickness cohesive interface elements across different phases. Further details on the CZM approach and the insertion of cohesive interface elements can be found in [PERSON] et al. [52]. Similarly, for the phase-field fracture simulations, surface continuity is maintained across material phases, although in this case it is not strictly necessary. As the interface thickness is set to 0.6 mm, meshing the model with an average element length of 0.05 mm ensures that there are sufficient elements within the interface zone to accurately capture the gradient terms of the phase field damage. As shown in Figure 4, all 2D specimens are subjected to a displacement of 0.05 mm in the \(x\) direction. Figure 4: **(a)** Benchmark model illustrating three distinct concrete mesoscale microstructure phases. **(b)** Applied boundary conditions used in all simulations conducted in this study. **(c)** Finite element discretization employed for phase-field fracture simulations. **(d)** Finite element discretization applied for cohesive zone simulations. Material properties for the CZM and PFM simulations are given in Table 2 and Table 3, respectively. Inclusions, known for their significantly higher mechanical strength, typically prevent crack propagation, a common observation in many engineering materials. As a result, the damage characteristics of inclusions are not explicitly defined in the CZM simulations. To reflect this behaviour in PFM simulations, the failure strength of inclusions is deliberately set much higher than that of the other material phases. In this study, a staggered scheme is used for phase-field fracture simulations, where the phase-field and displacement fields are solved sequentially. The staggered scheme proposed by [PERSON] and [PERSON] [59] is used to solve the system of nonlinear equations for PFM simulations. In addition, the CPFM approach is integrated into the ABAQUS UEL subroutine using the open source ABAQUS codes from [PERSON] and [PERSON] [59], with modifications applied to the UEL subroutine to convert the code from the AT2 standard phase-field fracture model to the CPFM fracture model. These modifications are based on the equations presented in Table 1.  
Similar to other studies, such as [57; 60], the Abaqus Explicit solver is used for the CZM simulations. The reason for this choice is that the insertion of CIE between each two solid elements introduces a high degree of nonlinearity into the model, making it very difficult to use static, quasi-static, or implicit solvers. However, in this study, to ensure that dynamic effects are not present in the CZM simulation, two strategies are implemented: first, the displacement is applied gradually, and second, the time increment is kept very small (in this study is \(2.99\times 10^{-8}\)). With these strategies, the kinetic energy in the CZM simulations remains around 2% of the total energy of CZM simulation. Moreover, reaction force-displacement curves are derived by averaging the reaction forces of the nodes along the left vertical edge, and stress-strain curves are obtained by dividing the reaction force by the element edge length of the benchmark model.