The simulation of spherical microphone arrays is commonly performed in the frequency domain using spherical harmonic representation of the spatial transfer functions. Each modal spectrum is described by the spherical Hankel function and its derivative. Although the resulting simulation is accurate in the frequency domain, the corresponding time-domain signal obtained by the inverse discrete Fourier transform (IDFT) exhibits temporal aliasing and smearing. Moreover, evaluating the spherical Hankel functions at a larger number of frequencies can be a computational bottleneck. In this paper, we propose a timedomain approach, where each modal transfer function is implemented as a parallel combination of IIR filters and a single FIR filter. The poles of the IIR filters correspond to the roots of the spherical Hankel functions' derivatives, which can be pre-computed for a given radius. The FIR filter coefficients are obtained by the least-squares solution, where the squared spectrum errors are minimized at control frequencies. While the number of poles are fixed for each harmonic order, the FIR length is a free design parameter, with which the simulation accuracy can be adjusted. The presented approach is compared numerically with previously proposed time-domain methods and the frequency-domain modeling.