Toughening strategies are needed to design high performance secondary bonded joints. Although adhesively bonded joints present a promising solution for cost-saving assembly, they are however among the most vulnerable parts in composite structures. Due to the inherent heterogeneity of composite substrates and the infeasibility of ideal surface pretreatments, the bonding properties of composite joints feature a spatial variability that can result in a wide variety of debonding phenomenology, such as crack tip transfer, crack bifurcation or ligament bridging. In the present work, we apply a cohesive zone model with spatially stochastic interfacial properties to model the physical behaviors that take place during debonding of an adhesively bonded composite joint. The joint distribution of the two essential bonding properties, namely, the critical separating stress σm and the strain energy release rate Gc, follows a bivariate lognormal probability distribution. The proposed model captures the mechanisms of crack-tip transfer and the failure of bridging adhesive ligaments between the composite adherends, as well as the associated toughening of joints. We find that the extrinsic dissipation attributed to extra debonding areas and plastic-damage mechanisms of ligaments may be one of the primary causes for the thickness effect of bondline on the fracture toughness, commonly observed for adhesively bonded joints. Through a parametric study, we discuss the effect of variability in interfacial properties on the enhancement of toughness of bonded composite joints. The outcomes of this work may provide insights into the fundamental failure mechanisms of secondary bonded composite joints.