The production of Z boson pairs in proton–proton ( (pp) ) collisions, (pp\to \left(Z/{\gamma }^{\ast }\right)\left(Z/{\gamma }^{\ast }\right)\to 2\ell 2{\ell }^{\prime }), where (\ell,{\ell }^{\prime }=e) or (\mu ), is studied at a center-of-mass energy of 13 (\phantom{\rule{0.166667em}{0ex}}\text{TeV}) with the CMS detector at the CERN LHC. The data sample corresponds to an integrated luminosity of 137 (\phantom{\rule{0.166667em}{0ex}}{\text{fb}}^{-1}), collected during 2016–2018. The (ZZ) production cross section, ({\sigma }{\text{tot}}\left(pp\to ZZ\right)=17.4±0.3\phantom{\rule{0.166667em}{0ex}}\text{(stat)}±0.5\phantom{\rule{0.166667em}{0ex}}\text{(syst)}±0.4\phantom{\rule{0.166667em}{0ex}}\text{(theo)}±0.3\phantom{\rule{0.166667em}{0ex}}\text{(lumi) pb}), measured for events with two pairs of opposite-sign, same-flavor leptons produced in the mass region (60<{m}{{\ell }^{+}{\ell }^{-}}<120\phantom{\rule{0.166667em}{0ex}}\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 (ZZZ) and (ZZ\gamma ) couplings.