The standard finite element implementation of intrinsic cohesive zone models (CZMs) based on the penalty method exhibits a distinct lack of numerical stability and/or convergence for stiff cohesive laws. This lack of stability is typically observed in the form of spurious oscillations in the normal and tangential tractions recovered at the cohesive interface. In this paper, we will present a robust, stabilized finite element formulation for CZMs that remedies traction oscillations, thus ensuring stability and convergence for any value of initial cohesive stiffness. A key advantage of the proposed formulation is that it generalizes the Nitsche's method for modeling cohesive fracture with a large initial cohesive stiffness, thus enabling the implementation of intrinsic and extrinsic CZMs in a unified and variationally consistent manner. We present several numerical examples to demonstrate the stability, convergence and accuracy of the proposed formulation in two-dimensions. First, we will verify the accuracy using simple patch tests considering uniaxial tension, compression and shear loadings. Second, we will demonstrate the lack of spurious traction oscillations at cohesive interfaces of rectangular beams loaded under shear and three-point bending. To demonstrate the stability issues related with the spurious traction oscillation, we consider both isotropic as well as anisotropic CZMs, wherein the normal and tangential cohesive stiffness values are different. Our numerical results for high stiffness cases clearly show that the proposed formulation yields a smooth oscillation-free traction profile and ensures stability, whereas the standard formulation suffers from instability and/or convergence issues. © 2018 by DEStech Publications, Inc. All rights reserved.