The production of Z boson pairs in proton-proton (pp) collisions, pp $\to$ (Z/$\gamma^*$)(Z/$\gamma^*$) $\to$ 2$\ell$ 2$\ell'$, where $\ell,\ell'$ = e or $\mu$, is studied at a center-of-mass energy of 13 TeV with the CMS detector at the CERN LHC. The data sample corresponds to an integrated luminosity of 137 fb$^{-1}$, collected during 2016-2018. The ZZ production cross section, $\sigma_{\text{tot}}$ (pp $\to$ ZZ) = 17.2 $\pm$ 0.3 (stat) $\pm$ 0.5 (syst) $\pm$ 0.4 (theo) $\pm$ 0.3 (lumi) pb, measured for events with two pairs of opposite-sign, same-flavor leptons produced in the mass region 60 $\lt m_{\ell^+\ell^-} \lt $ 120 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.