Statistically simulated time series of wave parameters are required for many coastal and offshore engineering applications, often at the resolution of approximately 1 h. Various studies have relied on autoregressive moving-average (ARMA) processes to simulate synthetic series of wave parameters in a Monte Carlo sense. However, accurately representing inter-series dependencies has remained a challenge. In particular, the relationship between wave height and period statistics is complex, due to the limiting steepness condition. Here, we present a new simulation method for joint time series of significant wave height, mean zero-crossing periods and a directional regime variable. The latter distinguishes between northern and southwestern waves. The method rests on several model components which include renewal processes, Fourier series with random coefficients, ARMA processes, copulas and regime-switching. A particular feature is a data-driven estimate for a wave height-dependent limiting wave steepness condition which is used to facilitate copula-based dependence modeling. The method was developed for and applied to a data set in the Southern North Sea. For this site, the method could simulate time series with realistic annual cycles and inter-annual variability. In the time series data, the bivariate distribution of significant wave height and mean zero-crossing period was well represented. An influence of the directional regime on the bivariate distribution could also be modeled. However, the influence was not as strong in simulated data as in observed data. Finally, simulated series captured duration and inter-arrival time of storm events well. Potential applications for output of the simulation method range from the assessment of risks to the planning and budgeting of operations in coastal and offshore environments.