This paper reports on a search for an extended scalar sector of the standard model, where a new CP-even (odd) boson decays to a Z boson and a lighter CP-odd (even) boson, and the latter further decays to a b quark pair. The Z boson is reconstructed via its decays to electron or muon pairs. The analysed data were recorded in proton-proton collisions at a center-of-mass energy $\sqrt{s} = $ 13 TeV, collected by the CMS experiment at the LHC during 2016, corresponding to an integrated luminosity of 35.9 fb$^{-1}$. Data and predictions from the standard model are in agreement within the uncertainties. Upper limits at 95% confidence level are set on the production cross section times branching fraction, with masses of the new bosons up to 1000 GeV. The results are interpreted in the context of the two-Higgs-doublet model.