Under varying climate conditions, cracks are commonly observed in bolted joints, owing to the shrinkage of wood and confinement from slotted-in steel plates and bolts. A three-dimensional finite element model was developed to investigate the mechanical behavior of bolted glulam joints with initial cracks. Wood foundation was prescribed in the model to simulate the local crushing behavior of wood surrounding the bolts. The behavior of wood in compression and the foundation were defined as transversely isotropic plastic in the software package ANSYS. Cohesive zone model was applied in the numerical analysis to consider the propagation of initial cracks and brittle failure of wood in the bolted joints under tension load. The numerical model was validated by the experiments conducted on full-scale specimens and it is indicated that the numerical model has good ability in predicting the failure modes and capacity of tension joints with local cracks. To further investigate the influence of crack number, length and locations, a parametric study was conducted with the verified model. Moreover, to study the effects of cracks on the behavior of bolted joints with different failure modes, another bolted joint including bolts with different strength grades and diameters was designed and analyzed in the parametric study, which was expected to have bolt yielding failure mode. It was found that the initial cracks can decrease the capacity and initial stiffness of tension joints by up to 16.5 and 34.8%, respectively.