Modelling multiple hazard interrelations remains a challenge for practitioners. This article primarily focuses on the interrelations between pairs of hazards. The efficacy of six distinct bivariate extreme models is evaluated through their fitting capabilities to 60 synthetic datasets. The properties of the synthetic datasets (marginal distributions, tail dependence structure) are chosen to match bivariate time series of environmental variables. The six models are copulas (one non-parametric, one semi-parametric, four parametric). We build 60 distinct synthetic datasets based on different parameters of log-normal margins and two different copulas. The systematic framework developed contrasts the model strengths (model flexibility) and weaknesses (poorer fits to the data). We find that no one model fits our synthetic data for all parameters but rather a range of models depending on the characteristics of the data. To highlight the benefits of the systematic modelling framework developed, we consider the following environmental data: (i) daily precipitation and maximum wind gusts for 1971 to 2018 in London, UK, and (ii) daily mean temperature and wildfire numbers for 1980 to 2005 in Porto District, Portugal. In both cases there is good agreement in the estimation of bivariate return periods between models selected from the systematic framework developed in this study. Within this framework, we have explored a way to model multi-hazard events and identify the most efficient models for a given set of synthetic data and hazard sets.
