For gauged watersheds, flood frequency analysis has been commonly done in one of three ways: (i) univariate frequency analysis of observed peak discharges at a gauged site; (ii) at-a-site bivariate frequency analysis of correlated flood variables, i.e., through the dependence of peak discharge and flood volume; flood volume and flood duration with the assumption of independence between peak discharge and flood duration; (iii) at-a-site trivariate flood frequency analysis of flood variables with the use of conventional trivariate distributions, symmetric 3-dimensional Archimedean copula, and asymmetric 3-dimensional Archimedean copula. However, for the three cases discussed above, there are some limitations as, for example: (i) the univariate analysis may not appropriately represent the flood risk; (ii) it may not be valid to assume the independence between peak discharge and flood duration; (iii) the chosen biviariate Archimedean copula may not be extended to higher-dimensional Archimedean copula (d >= 3), and (iv) there is parametric limitation for asymmetric 3-dimensional Archimedean copula. To address the aforementioned limitations, assuming univariate flood variables being Independently, Identically Distributed (I.I.D.) random variables, this paper presents a procedure to investigate the at-site dependence of flood variables (i.e., peak discharge, flood volume, and flood duration) based on the vine-copula theory. The proposed procedure is verified using experimental and natural watersheds in the U.S.