Not yet reviewed by Pith; the record is open.

Hundreds of millions of supermassive black hole binaries are expected to contribute to the gravitational-wave signal in the nanohertz frequency band. Their signal is often approximated either as an isotropic Gaussian stochastic background with a power-law spectrum, or as an individual source corresponding to the brightest binary. In reality, the signal is best described as a combination of a stochastic background and a few of the brightest binaries modeled individually. We present a method that uses this approach to efficiently create realistic pulsar timing array datasets using synthetic catalogs of binaries based on the Illustris cosmological hydrodynamic simulation. We explore three different properties of such realistic backgrounds which could help distinguish them from those formed in the early universe: i) their characteristic strain spectrum; ii) their statistical isotropy; and iii) the variance of their spatial correlations. We also investigate how the presence of confusion noise from a stochastic background affects detection prospects of individual binaries. We calculate signal-to-noise ratios of the brightest binaries in different realizations for a simulated pulsar timing array based on the NANOGrav 12.5-year dataset extended to a time span of 15 years. We find that $\sim$6% of the realizations produce systems with signal-to-noise ratios larger than 5, suggesting that individual systems might soon be detected (the fraction increases to $\sim$41% at 20 years). These can be taken as a pessimistic prediction for the upcoming NANOGrav 15-year dataset, since it does not include the effect of potentially improved timing solutions and newly added pulsars.