SummaryPair-copula construction methodology has been explored to model the dependence structure between net storm event depth (R), maximum wet periods' depth (M), and the total wet periods' duration (L), noting that the total storm event depth is RT = R + M. Random variable R was used instead of RT in order to avoid physical boundary effects due to the condition of RT ⩾ M. The flexibility of pair-copula construction allowed the examination of 11 bivariate copulas at the three bivariate stages of the three-dimensional (3D) copula. For 21 years of hourly rainfall data from Cook County, Illinois, USA, examined, three different copulas were found suitable for the bivariate stages. For the internal storm event structure, a Geometric distribution was used to model the net event duration, defined as the difference between the total duration (D) and L. A two-parameter Poisson model was adopted for modelling the distribution of the L wet periods within D, and the first-order autoregressive Lognormal model was applied for the distribution of RT over the L wet periods. Incorporation of an inter-event (I) sub-model completed the continuous rainfall simulation scheme. The strong seasonality in the marginal and dependence model parameters was captured using first harmonic Fourier series, thus, reducing the number of parameters. Polynomial functions were fitted to the internal storm event model parameters which did not exhibit seasonal variability. Four hundred simulation runs were carried out in order to verify the developed model. Kolmogorov-Smirnov (KS) tests found the hypothesis that the observed and simulated storm event quantiles come from the same distribution cannot be rejected at the 5% significance level in nearly all cases. Gross statistics (dry probability, mean, variance, skewness, autocorrelations, and the intensity-duration-frequency (IDF) curves) of the continuous rainfall time series at several aggregation levels were very well preserved by the developed model.