Reliability analysis of large coupled systems is a challenging task in dam safety. Not only the number of random variables is increased, but also a complex numerical model makes the simulations computationally expensive. To overcome this problem, an efficient analytical model should be developed which does not affect the accuracy of reliability analysis. In this paper, a gravity dam-foundation-reservoir coupled system is taken as a vehicle for the numerical simulations. Epistemic (i.e. material) uncertainties, as well as the ground motion record-to-record variability (i.e. aleatory) are incorporated in the simulations. A new reliability technique is proposed and its efficiency is compared to the Latin Hypercube Sampling (LHS). In this technique, a continuous probability density function is replaced by a finite number of fractional moments. Furthermore, an improved bivariate dimension reduction method is developed for dynamic systems. Several parametric simulations are performed to investigate the impact of the reservoir water level and the dam size. Results of the proposed method are in good agreement with the LHS. Finally, a response surface meta-model is proposed to estimate the dam behavior as a function of ground motion intensity measures.