We produce 1000 realizations of synthetic clustering catalogues for each type of the tracers used for the baryon acoustic oscillation and redshift space distortion analysis of the Sloan Digital Sky Surveys-IV extended Baryon Oscillation Spectroscopic Survey final data release (eBOSS DR16), covering the redshift range from 0.6 to 2.2, to provide reliable estimates of covariance matrices and test the robustness of the analysis pipeline with respect to observational systematics. By extending the Zel’dovich approximation density field with an effective tracer bias model calibrated with the clustering measurements from the observational data, we accurately reproduce the two- and three-point clustering statistics of the eBOSS DR16 tracers, including their cross-correlations in redshift space with very low computational costs. In addition, we include the gravitational evolution of structures and sample selection biases at different redshifts, as well as various photometric and spectroscopic systematic effects. The agreements on the auto-clustering statistics between the data and mocks are generally within 1σ variances inferred from the mocks, for scales down to a few h−1Mpc in configuration space, and up to 0.3hMpc−1 in Fourier space. For the cross correlations between different tracers, the same level of consistency presents in configuration space, while there are only discrepancies in Fourier space for scales above 0.15hMpc−1. The accurate reproduction of the data clustering statistics permits reliable covariances for multi-tracer analysis.