RAID systems are ubiquitously deployed in storage environments, both as standalone storage solutions and as fundamental components of virtualized storage platforms. Accurate models of their performance are crucial to delivering storage infrastructures that meet given quality of service requirements. To this end, this paper presents a flexible fork-join queueing simulation model of RAID systems that are composed of zoned disk drives and which operate under RAID levels 01 or 5. The simulator takes as input I/O workloads that are heterogeneous in terms of request size and that exhibit burstiness, and its primary output metric is I/O request response time distribution. We also study the effects of heavy workload, taking into account the request-reordering optimizations employed by modern disk drives. All simulation results are validated against device measurements and compared with existing analytical queueing network models for the development of the models.