Integrated safety analysis of hazardous process facilities calls for an understanding of both stochastic and topological dependencies, going beyond traditional Bayesian Network (BN) analysis to study cause-effect relationships among major risk factors. This paper presents a novel model based on the Copula Bayesian Network (CBN) for multivariate safety analysis of process systems. The innovation of the proposed CBN model is in integrating the advantage of copula functions in modelling complex dependence structures with the cause-effect relationship reasoning of process variables using BNs. This offers a great flexibility in probabilistic analysis of individual risk factors while considering their uncertainty and stochastic dependence. Methods based on maximum likelihood evaluation and information theory are presented to learn the structure of CBN models. The superior performance of the CBN model and its advantages compared to traditional BN models are demonstrated by application to an offshore managed pressure drilling case study. (C) 2016 Elsevier Ltd. All rights reserved.