This paper documents the sixteenth data release (DR16) from the Sloan Digital Sky Surveys; the fourth and penultimate from the fourth phase (SDSS-IV). This is the first release of data from the southern hemisphere survey of the Apache Point Observatory Galactic Evolution Experiment 2 (APOGEE-2); new data from APOGEE-2 North are also included. DR16 is also notable as the final data release for the main cosmological program of the Extended Baryon Oscillation Spectroscopic Survey (eBOSS), and all raw and reduced spectra from that project are released here. DR16 also includes all the data from the Time Domain Spectroscopic Survey (TDSS) and new data from the SPectroscopic IDentification of ERosita Survey (SPIDERS) programs, both of which were co-observed on eBOSS plates. DR16 has no new data from the Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey (or the MaNGA Stellar Library "MaStar"). We also preview future SDSS-V operations (due to start in 2020), and summarize plans for the final SDSS-IV data release (DR17).

The Baryon Oscillation Spectroscopic Survey of SDSS-III has provided an unprecedentedly large sample of Damped \lya systems (DLAs), the largest repositories of neutral hydrogen in the Universe. This DLA sample has been used to determine the DLA bias factor from their cross-correlation with the \lya forest absorption in \cite{FontRibera2012,Perez2018}, showing that DLAs are associated with relatively massive halos. However, the low resolution and signal-to-noise of BOSS spectra do not allow precise measurements of the DLA metal lines. We define a metal strength parameter, $S$, based on combining equivalent widths of 17 metal lines, that can be measured with an optimal signal-to-noise ratio for individual DLAs in BOSS spectra, allowing for the classification of the DLA population into subgroups of different $S$. We present the distribution of this DLA metal strength and the dependence of its mean value on $N_{\rm HI}$ and redshift. We search for systematic effects and variations in the catalogue purity by examining the dependence of the $S$ distribution on the spectral signal-to-noise and the estimated error on $S$. A catalogue of DLAs with measured equivalent widths for the selected 17 metal lines and the value of $S$ is made publicly available, which will be used to measure the dependence of the DLA bias factor on the $S$ parameter. The relation of the metal strength on the gas metal abundances and velocity dispersion can be constrained by studying the stacked metal absorption spectra of DLAs as a function of $S$, allowing for future determinations of the dependence of the bias factor on the metallicity and velocity dispersion of DLAs.