The mechanical behaviour and failure properties of adhesive joints are influenced not only by adhesive properties but also by interface characteristics. Various modelling methods have been proposed to predict the failure behaviour of adhesive joints. Among those methods, cohesive zone modelling, which consists of introducing a cohesive law in the numerical models to reproduce the adhesive layer, is widely used to simulate the behaviour of adhesive joints. The cohesive parameters of the traction–separation laws in each pure mode should be determined for application of this technique. It is known that the cohesive parameters are affected by both adhesive properties and interface characteristics so, in most cases, these parameters are iteratively obtained by matching numerical results to experimental results. In this paper, a systematic procedure for the determination of the cohesive parameters is proposed by introducing an optimization technique (design of experiment and the kriging metamodel). The first step is to identify design variables affecting a response of a system. In this study, the cohesive parameters are selected as the design variables. Once the design variables are identified, design ranges or sampling ranges of the design variables are determined in order to result in a near optimal response with the determined sampling ranges. If the design variables and the sampling ranges are obtained then the design of experiment is performed to determine appropriate sampling points in the sampling ranges and a minimum number of simulations. Then, the load difference between numerical and experimental load–displacement results at several values of displacements are defined as an error, and the kriging metamodel, which surrogates the error, is constructed based on the sampling points. Cohesive parameters are determined by applying the numerical nonlinear optimization algorithm to minimize the error. The proposed procedure is applied to a co-cured Single Leg Bending (SLB) joint under mode I and mode II-dominant modes and the cohesive parameters are determined. From these parameters, the mixed-mode cohesive zone model is constructed and applied to the simulation of the co-cured SLB joint under varying mode-mixities. Simulation results are compared with the experimental results. It is shown that the failure behaviour of the SLB joint is well-described by using the mixed-mode cohesive zone model with the determined parameters.