Uncertainty exists widely in engineering practice. An engineering system may have multiple failure criteria. In the current paper, system reliability analysis with multiple failure modes under both epistemic and aleatory uncertainties is presented. Epistemic uncertainty is modelled using p-boxes, while aleatory uncertainty is modelled using probability distributions. A first-order reliability method is developed and non-linear performance functions are linearized by the sampling method instead of the commonly used Taylor’s expansion at the most probable point. Furthermore, multiple failure modes in a system are often correlated because they depend on the same uncertain variables. In order to consider these correlated failure modes, the methods proposed by Feng and Frank are extended in this paper in order to calculate the joint probability of failure for two arbitrary failure modes under both aleatory and epistemic uncertainties. The Pearson correlation coefficient of two arbitrary failure modes is determined by the sampling method. Since two types of uncertainty exist in the system, the probability of system failure is an interval rather than a point value. The probability of failure of the system can be obtained by the combination of the extension ‘narrow’ bound method and the interval arithmetic. A numerical example is presented to demonstrate the applicability of the proposed method.