We study the critical behavior of the three-dimensional (3D) Gross-Neveu (GN) model with $N_f$ Dirac fermionic flavors and quartic interactions, at the chiral ${\mathbb Z}_2$ transition in the massless ${\mathbb Z}_2$-symmetric limit. For this purpose, we consider a lattice GN model with staggered Kogut-Susskind fermions and a scalar field coupled to the scalar bilinear fermionic operator, which effectively realizes the attractive four-fermion interaction. We perform Monte Carlo (MC) simulations for $N_f=4,8,12,16$. By means of finite-size scaling analyses of the numerical data, we obtain estimates of the critical exponents that are compared with the large-$N_f$ predictions obtained using the continuum GN field theory. We observe a substant...