The production of a Z boson, decaying to two charged leptons, in association with jets in proton-proton collisions at a centre-of-mass energy of 13 TeV is measured. Data recorded with the CMS detector at the LHC are used that correspond to an integrated luminosity of 2.19 fb$^{-1}$. The cross section is measured as a function of the jet multiplicity and its dependence on the transverse momentum of the Z boson, the jet kinematic variables (transverse momentum and rapidity), the scalar sum of the jet momenta, which quantifies the hadronic activity, and the balance in transverse momentum between the reconstructed jet recoil and the Z boson. The measurements are compared with predictions from four different calculations. The first two merge matrix elements with different parton multiplicities in the final state and parton showering, one of which includes one-loop corrections. The third is a fixed-order calculation with next-to-next-to-leading order accuracy for the process with a Z boson and one parton in the final state. The fourth combines the fully differential next-to-next-to-leading order calculation with next-to-next-to-leading logarithm resummation and parton showering.