Clustered bivariate survival data arise in various fields, such as biology and medicine, when individuals in a dataset are clustered and exhibit two survival outcomes. Recently, the joint frailty-copula model was proposed to analyze clustered bivariate survival outcomes by accommodating the between-cluster heterogeneity via a shared frailty term. In this model, researchers fitted the baseline hazard functions via the nonparametric model, the spline model, or the Weibull model. However, when a population has extremely large survival time, the baseline hazard functions are better modeled by a heavy-tailed distribution. In this paper, we adopt the Pareto type I distribution for the joint frailty-copula model, which is one of the most popular heavy-tailed distributions. We show that the moments of the Pareto type I joint frailty copula model diverge to infinity owing to the heavy right-tail. We develop statistical inference methods based on three types of censoring schemes: (i) bivariate random censoring, (ii) semi-competing risks, and (iii) competing risks. We develop maximum likelihood estimation procedures, and make our computational tools available for users. Simulations are performed to check the accuracy of the proposed method. We finally analyze a real dataset for illustration.