This paper presents a measurement of the production cross-section of a Z boson in association with b- or c-jets, in proton–proton collisions at $\sqrt{s} = 13$ TeV with the ATLAS experiment at the Large Hadron Collider using data corresponding to an integrated luminosity of 140 $\hbox {fb}^{-1}$. Inclusive and differential cross-sections are measured for events containing a Z boson decaying into electrons or muons and produced in association with at least one b-jet, at least one c-jet, or at least two b-jets with transverse momentum $p_\text {T} > 20$ GeV and rapidity $|y| < 2.5$. Predictions from several Monte Carlo generators based on next-to-leading-order matrix elements interfaced with a parton-shower simulation, with different choices of flavour schemes for initial-state partons, are compared with the measured cross-sections. The results are also compared with novel predictions, based on infrared and collinear safe jet flavour dressing algorithms. Selected $Z + \ge 1~c$-jet observables, optimized for sensitivity to intrinsic-charm, are compared with benchmark models with different intrinsic-charm fractions.