The production of Z boson pairs in proton–proton (${\mathrm{p}} {\mathrm{p}} $) collisions, ${{\mathrm{p}} {\mathrm{p}} \rightarrow ({\mathrm{Z}}/\gamma ^*)({\mathrm{Z}}/\gamma ^*) \rightarrow 2\ell 2\ell '}$, where ${\ell ,\ell ' = {\mathrm{e}}}$ or ${{\upmu }}$, is studied at a center-of-mass energy of 13$\,\text {TeV}$ with the CMS detector at the CERN LHC. The data sample corresponds to an integrated luminosity of 137$\,\text {fb}^{-1}$, collected during 2016–2018. The ${\mathrm{Z}} {\mathrm{Z}} $ production cross section, $\sigma _{\text {tot}} ({\mathrm{p}} {\mathrm{p}} \rightarrow {\mathrm{Z}} {\mathrm{Z}} ) = 17.4 \pm 0.3 \,\text {(stat)} \pm 0.5 \,\text {(syst)} \pm 0.4 \,\text {(theo)} \pm 0.3 \,\text {(lumi)} \text { pb} $, measured for events with two pairs of opposite-sign, same-flavor leptons produced in the mass region ${60< m_{\ell ^+\ell ^-} < 120\,\text {GeV}}$ is consistent with standard model predictions. Differential cross sections are also measured and agree with theoretical predictions. The invariant mass distribution of the four-lepton system is used to set limits on anomalous ${\mathrm{Z}} {\mathrm{Z}} {\mathrm{Z}} $ and ${{\mathrm{Z}} {\mathrm{Z}} \gamma }$ couplings.