A three-dimensional model of polydisperse reactive sedimentation is developed by means of a multilayer shallow water approach. The model consists of a variety of solid particles of different sizes and densities, and substrates diluted in water, which produce biochemical reactions while the sedimentation process occurs. Based on the Masliyah–Lockett–Bassoon settling velocity, compressibility of the sediment and viscosity of the mixture, the system of governing equations is composed by non-homogeneous transport equations, coupled to a momentum equation describing the mass-average velocity. Besides, the free-surface depicted by the total height of the fluid column is incorporated and fully determined through the multilayer approach. A finite volume numerical scheme on Cartesian grids is proposed to approximate the model equations. Numerical simulations of the denitrification process exemplify the performance of the numerical scheme and model under different scenarios and bottom topographies.