  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 $\sigma$ variances
inferred from the mocks, for scales down to a few $h^{-1}\,{\rm Mpc}$ in
configuration space, and up to $0.3\,h\,{\rm Mpc}^{-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.15\,h\,{\rm Mpc}^{-1}$. The accurate reproduction of
the data clustering statistics permits reliable covariances for multi-tracer
analysis.
